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Abstract 

This  report  treats  the  numerical  computation  of  cumulative  distributions  of 
random  variables  occurring  primarily  in  communications,  radar,  and  optics  when 
their  moment-generating  or  probability-generating  functions  are  known.  The 
cumulative  distribution  of  a  continuous  random  variable  is  expressed  as  a  Laplace 
inversion  integral  of  its  moment-generating  function,  that  of  an  integer-valued 
random  variable  as  a  contour  integral  that  arises  from  Cauchy's  theorem  and  whose 
integrand  involves  the  probability-generating  function.  These  integrals  are 
evaluated  by  numerical  quadrature  along  contours  in  the  complex  plane  chosen  for 
efficiency  and  convenience.  Applications  include  radar  detection  probabil ities 
with  fading  and  unfading  signals  and  fixed-threshold  and  constant-false-alarm- 
rate  receivers;  the  distributions  of  the  integrated  output  of  a  linear  rectifier 
and  of  the  filtered  output  of  a  quadratic  rectifier;  the  error  probability  in  a 
binary  symmetric  communication  channel  with  intersymbol  and  cochannel 
interference;  the  distribution  of  shot  noise;  the  distributions  of  the  numbers  of 
electrons  emerging  from  photoelectric  detectors,  photomultipliers,  and  avalanche 
diodes;  and  significance  probabilities  in  statistical  rank  tests. 
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I.  Distributions  of  Continuous  Random  Variables 


1.  Computing  Cumulative  Distributions  from  Moments  and  Moment -Generating 
Functions 

A  receiver  that  is  to  detect  a  particular  type  of  signal,  as  in  radar  and 
communications,  processes  its  input  in  such  a  way  as  to  generate  a  random 
variable  that  we  shall  denote  by  v.  If  v  exceeds  a  decision  level  V,  the 
receiver  decides  that  the  signal  is  present.  The  efficacy  of  the  receiver  is 
measured  by  its  probability  of  detection,  which  is  the  probability 

q+(V)  =  Pr  (v  >  V)  =  J  p(v)  dv  0.1) 

V 


that  it  makes  an  affirmative  decision  when  its  input  consists  of  the  signal  in 
the  presence  of  random  noise  of  a  specified  character.  Here  p(v)  is  the 
probability  density  function  (p.d.f.)  of  the  random  variable  v,  and  q+(V)  is  the 
complementary  cumulative  distribution  function  (c.c.d.f.).  The  probability  q+(V) 
may  be  needed  as  a  function  of  the  input  signal-to-noise  ratio.  Also  relevant  is 
the  false-alarm  probability  that  v  >  V  when  only  noise  is  present. 

Depending  on  the  nature  of  the  signal  and  the  noise  and  on  the  structure  of 
the  receiver,  probabilities  such  as  q+(V)  may  be  very  difficult  to  calculate.  In 
many  situations,  however,  it  is  not  so  difficult  to  determine  the  moment¬ 
generating  function  (m.g.f.)  h(z)  of  the  random  variable  v,  which  is  the  Laplace 
transform  of  its  p.d.f.  p(v), 


h(z) 


(1  .2) 


In  other  problems,  although  a  closed  form  for  the  m.g.f.  h(z)  may  be  unknown,  one 
can  determine  the  moments 
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p(v)dv,  n  =  1 ,  2,  .... 


(1.3) 


of  the  random  variable  v,  which  are  the  coefficients  of  the  power  series 


h(z)  =  1  +  Pn(-z)n/n!  =  1  +  ^  an(-z)n 


(1.U) 


representing  the  m.g.f.  in  the  neighborhood  of  the  origin.  The  principal  topic 
of  this  first  part  of  our  report  is  the  determination  of  the  c.c.d.f.  q+(V)  from 
the  m.g.f.  h(z) . 

Detectors  of  light  signals,  as  in  optical  communications  and  radar,  often 
base  their  decisions  about  the  presence  or  absence  of  a  signal  on  the  number  k  of 
photoelectrons  counted  at  the  output  of  a  device  such  as  a  photoelectric  cell,  a 
photomultiplier,  or  an  avalanche  diode.  The  probability 


Qn  =  Pr  (k  >  n)  =  >  pR 


(1.5) 


that  this  number  k  exceeds  an  integer-valued  decision  level  n  is  useful  for 
characterizing  the  performance  of  the  detector.  Here  pk  is  the  probability  that 
exactly  k  electrons  are  counted.  Although  it  may  not  be  known  for  all  values  of 
k,  it  is  possible  in  many  problems  to  determine  the  probability  generating 
function  (p.g.f.)  h(z)  of  the  random  variable  k,  defined  as 


h(z)  =  >  pkz  . 


(1.6) 


Our  research  has  dealt  primarily  with  numerical  methods  for  calculating 
probabilities  such  as  those  in  (1.1)  and  (1.5)  from  the  generating  functions 
h(z).  Part  I  of  this  report  treats  the  distributions  of  continuous  random 
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variables  such  as  v;  Part  II  treats  the  distributions  of  integer-valued  random 
variables  such  as  k. 

(a)  The  Inversion  Integrals 

It  is  assumed  that  the  m.g.f.  h(z)  is  regular  at  least  in  a  vertical 
convergence  strip  parallel  to  and  containing  the  Im  z-axis  of  the  complex 
plane.  The  inverse  transform  to  (1.2)  is 


p(v)  = 


/c+i® 

.  h 


V7 

(z)e  dz/2-rri, 


(1.7) 


where  c  lies  in  the  convergence  strip  and  on  the  Re  z-axis.  The  cumulative 
distribution  function  (c.d.f.) 


q  (V)  =  Pr  (v  <  V)  =  I  p ( v )  dv 


can  be  obtained  by  integrating  (1.7),  provided  that  c  is  positive, 


(1.8) 


q"(V)  = 


-1  1/7 

z  h(z)e  dz/2iri, 


c  >  0. 


(1.9) 


Its  complement,  which  we  have  called  the  c. c.d.f.  or  the  probability  of 
detection ,  is  by  (1.1) 


q+(V)  =  1  -  q  (V)  =  -  /  z  1  h(z)eVz  dz/2Tri,  c  <  0.  (1.10) 


These  integrals  cannot  often  be  expressed  in  closed  form,  and  one  must  resort  to 
numerical  methods  for  evaluating  the  cumulative  distribution  and  its  complement. 

The  most  common  procedure  utilizes  the  Edgeworth  series.  One  defines  the 
cumulant  expansion  by 
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Tr. 


In  h(z)  =  ^  ]  t>n(-z)n  =  -V z  +  V2  o2z2  +  r(z), 


r(z)  =  /  v  <n^-z^  /n! > 


(1.11) 


where  V  =  E(v),  o  =  Var  v,  and  k  =  n!b  is  the  n-th  cumulant  of  v.  When  the 

n  n 


m.g.f.  h(z)  is  not  given  in  closed  form,  the  coefficients  b  =  k  /n!  of  this 

n  n 


series  can  be  obtained  from  the  coefficients  a  =  y  /n!  of  (1.4)  by  the  recurrent 

n  n 


relations 


"'n  +  l  an  +  1  n  +  1  ^  -1  rbrarn-1  -r  ' 


n  >  1  , 


(1.12) 


starting  with  b^  =  a,)  =  V. 


One  then  expands  the  m.g.f.  from  (1.11)  as 


h(z)  =  exp  (-Vz  +  *4  o2z2)  j  1  +  y  '  ^~kl^ 


(1.13) 


writes  out  the  terms  in  the  brackets  as  a  series  in  powers  of  (-z),  substitutes 


it  into  (1.7),  and  integrates  term  by  term  to  obtain  a  series  for  the  probability 


density  function  p(v).  This  series  is  in  turn  integrated  term  by  term  to 


produce,  after  appropriate  rearrangement,  the  well-known  Edgeworth  series  for  the 


c.d.f.  q  (V).  It  is  an  asymptotic  series  in  the  sense  that  although  its  terms  at 


first  decrease  in  magnitude,  the  series  diverges  when  carried  beyond  a  certain 


point.  Those  terms,  proportional  to  derivatives  of  the  error  function,  are 


oscillatory  functions  of  their  argument  (V-V)/o,  and  the  terms  of  high  order 


attain  large  positive  and  negative  values.  Computation  with  the  series  is 


difficult  to  manage,  therefore,  and  the  results  are  unreliable  in  the  tails  of 


■\ V  -W'/-'  ,v 
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the  distribution  where  |V  -  Vj  >>  a. 

(b)  The  Path  of  Steepest  Descent 

In  our  approach  one  avoids  the  expansion  in  (1.1 3).  which  is  the  cause  of 
all  the  trouble,  and  simply  integrates  (1.9)  or  (1.10)  numerically.  The 
numerical  quadrature  formulas  we  have  used  divide  the  integration  variable  into 
equal  intervals  or  steps  and  approximate  the  integral  as  a  weighted  sum  of  the 
values  of  the  integrand  at  the  ends  of  those  intervals.  Once  the  value  of  the 
integrand  has  decreased  below  a  certain  limit,  one  halts  the  summation. 

For  efficiency  one  would  like  to  minimize  the  number  of  such  steps  required 
for  a  given  accuracy,  and  one  therefore  desires  the  absolute  value  of  the 
integrand  to  drop  off  to  zero  as  rapidly  as  possible  along  the  contour  of 
integration.  The  contour  along  which  the  absolute  value  of  the  integrand 
decreases  most  rapidly  to  zero  is  called  its  path  of  steepest  descent  [1].  This 
path  crosses  the  Re  z-axis  at  a  saddlepoint  xQ  of  the  integrand,  which  is  a  root 
of  the  equation 

r(z)  =  =  ~~  In  h(z)  +  V  -  Z-1  =0,  d.u; 

dz  dz 

where 

T(z)  =  In  h(z)  +  Vz  -  In  (±z)  (1.15) 

is  the  "phase"  of  the  integrand  in  (1.9)  or  (1.10).  [(±z)  =  +z  for  q~(V)  in 
(1.9);  (±z)  =  -z  for  q+(V)  in  (1.10).]  Viewed  as  a  function  of  x  =  Re  z,  Y(z)  is 
minimum  at  the  saddlepoint  z  =  Xq.  Along  the  path  of  steepest  descent,  which  may 
have  several  branches,  Im  T(z)  is  constant.  On  the  principal  branch  crossing  the 
Re  z-axis  at  z  =  Xq,  Im  't'(z)  =  0.  In  many  problems  this  is  the  only  branch. 

Because  computing  the  path  of  steepest  descent  during  the  numerical 
integration  would  be  cumbersome,  one  prefers  to  approximate  it,  if  possible,  by  a 
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simple  curve  that  lies  close  enough  to  it  that  the  integrand  will  decrease  almost 


as  rapidly.  In  seeking  such  a  curve,  it  is  helpful  to  plot  paths  of  steepest 
descent  for  a  few  typical  values  of  the  parameters  of  the  distribution.  One  must 
first  determine  the  saddlepoint  Xq  where  the  path  crosses  the  Re  z-axis  by 
solving  (1.14).  One  can  often  do  so  most  expeditiously  by  Newton's  method,  which 
replaces  a  trial  value  of  Xq  at  each  stage  by 


vo  '  xo  ‘  >r(x0)’ 


Y'(x0) 


(1.16) 


primes  indicating  differentiation.  Because  of  the  convexity  U  of  the  phase  ¥(z) 
for  z  =  x  on  the  Re  z-axis,  this  method  converges  rapidly.  The  root  xQ  of  (1.14) 
is  sought  in  0  <  x  <  Y+  for  ( 1 . 9)  and  in  Y_  <  x  <  0  for  (1.10),  where  Y_  is  the 
rightmost  singularity  of  h(z)  on  the  negative  real  axis  and  Y+  is  the  leftmost 
singularity  of  h(z)  --  if  any  --  on  the  positive  real  axis.  If  it  is 
inconvenient  to  calculate  the  second  derivative  V"(z),  the  secant  method  can  be 
used  instead.  If  even  the  first  derivative  ¥'(z)  is  difficult  to  calculate,  one 
can  solve  the  equation 


Itn  Y(x0  +  ie)  =  0  (1.17) 

by  the  secant  method  for  a  sufficiently  small  positive  value  of  e,  or  one  can 
minimize  Y(z)  in  the  region  Y_  <  x  <  0  or  0  <  x  <  Y+  by,  for  instance,  the 
golden-section  algorithm  [2], 

One  traces  an  approximate  path  of  steepest  descent  by  starting  at  a  point  z1 
-  Xq  ♦  iS  a  small  distance  5  above  the  saddlepoint.  Given  the  k-th  point  zk  on 
the  approximate  path,  one  calculates  the  segment  Az  to  the  next  point  zk+1  = 
zk  +  Az  by  the  formula  [3] 

r*(z  ) 

z,  ,  -  z.  =  Az  =  -  6  r;:77 - rr,  (1.18) 

k+1  k  V'(z  ) 
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checking  the  value  of  Im  y(zk)  occasionally  to  see  whether  it  remains  close  to 
the  constant  value  C  for  the  path.  When  it  begins  to  deviate  unduly,  one  returns 
to  the  exact  path  by  solving  the  equation 

Im  ¥ (x+iy )  =  C 

by  Newton's  method  or  the  secant  method,  searching  in  x  when  the  slope 

|tan  (arg  Az)j  exceeds  1  and  in  y  when  it  is  less  than  1.  (On  the  main  branch, 

C  =  0.)  It  is  instructive  also  to  calculate  the  magnitude  exp  [Re  V(z)J  of  the 
integrand  to  see  whether  it  is  decreasing  with  satisfactory  rapidity  along  this 
path. 

If  the  main  branch  of  the  path  of  steepest  descent  heads  toward  infinity  in 
such  a  direction  that  the  originally  vertical  contour  in  (1.9)  or  (1.10)  can  be 
deformed  into  it  without  crossing  any  singularities  of  the  integrand,  even  at 
infinity,  all  is  well,  and  this  is  the  only  significant  branch.  In  some 
problems,  however,  the  path  along  which  Im  f(z)  =  0  may  become  parallel  to  the  Re 
z-axis  or  even  approach  it  asymptotically.  Such  a  behavior  indicates  that 
singularities  of  the  integrand  may  lie  at  infinity  and  that  the  path  of  steepest 
descent  possesses  several  or  even  an  infinite  number  of  branches.  These  come  in 
from  infinity,  turn  around,  and  pass  back  out  to  infinity.  On  each  such  branch 
lies  a  saddlepoint  of  the  integrand,  and  in  order  to  trace  it,  one  must  first 
locate  its  saddlepoint.  Again  one  solves  the  equation  V(z)  =  0  for  what  is  now 
a  complex  root  z.  Newton's  method  (1.16)  can  be  used,  with  Xq  replaced  by 
complex  trial  values  zQ,  but  the  selection  of  the  initial  trial  value  generally 
requires  preliminary  analysis  of  the  function  f(z)  in  the  complex  plane.  If  one 
starts  too  far  away,  Newton's  method  will  go  awry. 

The  path  of  steepest  descent  passes  through  each  saddlepoint  Zq  in  a 
direction  specified  by  that  of  the  complex  number  ^"(Zq).  By  starting  at  points 
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zQ  ±  i6y"(z0)/|y" ( zQ ) j 


at  a  small  distance  6  on  each  side  of  the  saddlepoint  Zq,  one  can  trace  that 
branch  of  the  path  of  steepest  descent  by  the  method  just  described.  Along  it, 

Im  'j'(z)  is  constant, 

Im  Vi z)  =  Im  V(zQ)  =  C, 
but  it  is  not  necessarily  zero. 

If  the  m.g.f.  h(z)  possesses  zeros  above  and  below  the  Re  z-axis,  the  path 
of  steepest  descent  through  the  saddlepoint  Xq  may  pass  through  all  or  some  of 
them  as  it  goes  off  to  infinity.  Between  each  pair  of  zeros  lies  another 
saddlepoint  of  the  integrand.  Even  more  complicated  configurations  of 
saddlepoints  and  path  segments  may  occur,  as  illustrated  by  Rice  [3],  and  finding 
a  suitable  contour  of  integration  may  be  most  difficult.  One  may  then  resort  to 
the  original  straight  vertical  contour  in  (1.9)  and  (1.10).  In  some  problems  the 
value  of  the  integrand  along  this  straight  path  may  oscillate  and  decrease  so 
slowly  that  other  methods  of  inverting  the  m.g.f.  must  be  sought.  Rice  [4]  has 
shown  how  in  certain  instances  one  can  make  the  integral  converge  more  rapidly  by 
subtracting  out  terms  whose  integrals  along  the  contour  are  known  exactly. 

It  may  happen  that  one  knows  the  moments  pn  of  the  random  variable  v,  but 
not  the  analytical  form  of  the  m.g.f.  h(z).  One  can  then  express  the  phase  Viz) 
as 

00 

'Viz)  =  bn(-z)n  +  Vz  -  In  (±z),  (1.19) 

n  =  1 

whose  coefficients  bn  =  <n/n!  can  be  determined  from  the  quantities  an  =  un/n!  by 
the  recurrence  in  (1.12).  One  stops  the  summation  in  (1.19)  when  the  term 
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t>n(-z)n  becomes  sufficiently  small.  Alternati vely  one  can  take 

H'(z)  =  In  h(z)  +  Vz  -  In  (±z) 


and  compute  the  m.g.f.  h(z)  from  the  power  series  in  (1.3).  In  many  problems  the 
magnitude  of  the  integrand  of  (1.9)  or  (1.10)  decreases  to  a  level  at  which  the 
integration  can  be  terminated  while  the  point  z  =  x  +  iy  still  lies  within  the 
domain  of  convergence  of  the  series  in  (1.3)  and  (1.19). 

(c)  Parabolic  Paths 

When,  as  often,  all  the  singularities  of  the  m.g.f.  lie  on  the  negative  Re 
z-axis,  the  path  of  steepest  descent  ultimately  goes  off  to  infinity  along  a  path 
resembling  a  parabola,  and  the  original  vertical  contour  in  (1.9)  or  (1.10)  can 
be  deformed  into  this  path  without  crossing  any  singularities  of  the  integrand. 
One  can  then  take  as  the  contour  of  integration  a  parabola  passing  through  the 
saddlepoint,  on  which 

z  =  xQ  +  J/2  <y^  +  iy  (1.20) 

witn  <  the  curvature  at  the  saddlepoint.  The  integrals  to  be  evaluated  can  be 
put  into  the  form 


q±(V) 


(l-i<y) 


dy/ir, 


(1  .21  ) 


with  T(z)  given  by  (1.15).  For  V  >>  E(v)  it  is  most  accurate  to  compute  q+(V), 
and  for  that  integration  Y_  <  xQ  <  0;  for  V  <<  E(v)  one  computes  q~(V),  for  which 
0  <  Xq  <  y+.  For  V  near  E(v)  it  does  not  much  matter  which  saddlepoint  one 
chooses . 

Designating  the  dominant  terms  in  the  phase  T(z)  far  from  the  saddlepoint  by 
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'j'(z),  one  calculates  the  curvature  <  of  the  approximating  parabola  at  y  =  0  by 
[5] 


(1.22) 


(1.23) 


The  curvature  <  is  not  needed  to  great  accuracy. 

The  numerical  integration  is  carried  out  by  the  trapezoidal  rule,  which 
Schwartz  [6]  and  Rice  [3]  have  shown  to  be  superior  for  infinite  integrals  of 
analytic  functions.  The  step  size  Ay  in  the  integration  variable  y  can  be  taken 
initially  as 

Ay  =  [Y<T(x0)]'^,  (1.24) 

which  with  Y  of  the  order  of  1  or  2  roughly  measures  the  width  of  the  integrand 
of  (1.21)  in  the  neighborhood  of  the  saddlepoint.  One  then  successively  halves 
the  step  size  Ay  until  the  values  of  the  integral  cease  changing  in  the  number  of 
figures  desired.  In  many  problems  the  number  of  significant  figures  roughly 
doubles  each  time  one  divides  Ay  by  2  [6].  The  numerical  integration  is 
conveniently  terminated  when  the  absolute  value  exp  [Re  ¥(2)]  of  the  last  term 
added  to  the  trapezoidal  sum  falls  below  cAy  times  the  absolute  value  of  the  sum 

itself,  with  e  some  suitably  small  number  depending  on  the  number  of  significant 

__  £ 

figures  desired.  We  found  e  =  10”°  adequate  for  most  purposes. 

(d)  Integration  Along  a  Path  of  Steepest  Descent 

When  the  parabola  is  not  a  sufficiently  close  approximation  to  the  path  of 
steepest  descent,  it  may  be  possible  to  find  some  other  simple  curve  along  which 
to  integrate;  an  example  will  be  given  in  Section  8.  If  the  path  of  steepest 
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descent  has  many  branches,  as  discussed  in  part  (b),  one  may  decide  to  integrate 
instead  along  a  close  approximation  to  each  branch  plotted  by  the  method 
described  there. 

That  approximate  path  consists  of  short,  straight  segments  Az  of  length  6 
connecting  a  sequence  of  points  zQ,  z1 ,  z2,  ...,  with  Az  given  by  (1.18).  One 
way  of  approximating  the  integral  along  any  given  semibranch  starting  at  zQ, 


Q  = 


dz/2ni , 


(1.25) 


is  to  replace  each  term  by  what  corresponds  to  the  trapezoidal  approximation, 


Re 


f  (z)dz/2iri 


f(z) 


=  Re  { L f ( )  +  f(zk+1)]  Az/Hiti}, 
=  exp  O(z)]. 


(1  .26) 


By  a  semibranch  we  mean  any  portion  of  the  path  beginning  at  a  saddlepoint  zQ  and 
extending  to  infinity;  the  branches  of  the  path  of  steepest  descent  above  the 
real  axis  consist  of  two  semibranches,  one  entering  zQ  from  infinity,  the  other 
leaving  it  in  the  opposite  direction.  The  contribution  of  the  entering  branch  is 
equal  to  the  negative  of  the  sum  in  (1.25). 

The  advantages  of  the  trapezoidal  rule  demonstrated  by  Schwartz  [6]  and  Rice 
[3]  for  integrals  of  analytic  functions  seem  not  to  carry  over  to  this 
situation.  Experience  has  shown  that  as  the  step  size  5  =  |Az|  is  halved,  the 
values  of  the  trapezoidal  sum  do  not  converge  with  the  rapidity  observed  when 
applying  the  trapezoidal  rule,  for  instance,  to  an  integral  such  as  (1.21).  We 
have  found  superior  convergence  when  instead  we  used  the  five-point  formula  of 
Birkhoff  and  Young  [73; 


1 1 


Re J  f(z)  dz/2iri  =  Re  1^7  {24f(zJ  ♦  U[f(z1#)  ♦  f(z,_J] 

z, 


60iri  - “k'  ''■*  v“k'  *'“k+1 

iAz .  ,  iAz. 


-  Cf(zk  -  *  f(zk  - 


(1  .27) 


where 


2k  •  (\  *  Vi)/2-  (,-28) 

The  error  Is  of  the  order  of  ; ’ ,  but  for  each  segment  of  the  contour  one  must 
compute  four  times  as  many  values  of  the  integrand  f(z)  =  exp[*F(z)]  as  for 
(1.26).  It  is  not  difficult  to  program  this  numerical  quadrature  in  connection 
with  the  procedure  embodied  in  ( 1 .18)  for  tracing  the  approximating  path  of 
steepest  descent. 

The  contribution  of  each  branch  of  the  path  of  steepest  descent  can  be 
roughly  estimated  by  the  saddlepoint  approximation  [3,  8] 

Q  =  Re  {[2it'F"(z0)]“^  exp  V(zQ) } ,  (1.29) 

where  zQ  is  the  saddlepoint  on  the  branch.  By  means  of  this  formula  one  can 
decide  whether  it  is  necessary  to  carry  out  the  numerical  integration  on  the 


branch. 


2.  Quadratic  Detection  of  Radar  Signals. 


(a)  Signals  of  Fixed  or  Random  Strength 

The  input  to  a  radar  receiver  is  passed  through  a  filter  matched  to  the 
shape  of  the  expected  echo  signal,  and  the  output  of  the  filter  is  applied  to  a 
quadratic  rectifier.  The  noise  in  the  input  is  white  and  Gaussian;  echoes  from 
a  target  may  also  be  present.  The  sampled  outputs  of  the  rectifier  in  a 
particular  range  bin  are  summed  in  M  successive  interpulse  intervals  to  produce  a 
detection  statistic  of  the  form 


M 


Z  = 


(2.1) 


J-1 

p 

where  r j  is  the  output  of  the  rectifier  in  the  j-th  interval.  When  suitably 
normalized,  Z  has  a  chi-squared  distribution  with  2M  degrees  of  freedom  under 
hypothesis  HQ  that  no  signal  is  present;  its  moment  generating  function  (m.g.f.) 
is 


hQ(u)  =  E(e”uZ 1  Hq)  =  (1  +  u)  M. 


(2.2) 


Under  hypothesis  H1  that  a  target  is  present,  an  echo  is  received  in  the 
j-th  interval  with  energy  and  signal-to-noise  ratio  =  2E-/ N,  1  £  j  £  Mf 

w  \J  y) 

where  N  is  the  unilateral  spectral  density  of  the  white-noise  background.  Then 
the  m.g.f.  of  Z  is 


h1  (u)  =  E(e  uZj  H.j )  =  (1  +  u)  M  exp  (-  . 

M 

s  =  y2  d2  «  y2  ^  dk2  =  et/n, 

k=1 


(2.3) 


where  Ej  is  the  total  received  energy.  The  probability  of  detection  is  given 
by  the  M-th  order  Q-function, 


8fl  -  Pr  (Z  >  Zj  8,) 


(2.«) 


where  Zq  is  the  decision  level  with  which  the  statistic  Z  is  compared  [9,  pp. 
215-219].  The  false-alarm  probability 


Q 


0 


(2.5) 


can  be  expressed  in  terms  of  the  incomplete  gamma  function.  A  saddlepoint 
approximation  was  applied  in  [8]  to  the  calculation  of  the  decision  level  Zq  and 
the  probability  Qd  of  detection  for  large  numbers  M  of  pulses  integrated;  typical 
results  are  tabulated  there  and  compared  with  the  exact  probabilities. 

When  the  radar  echoes  fade,  the  distribution  of  the  statistic  Z  must  be 
averaged  with  respect  to  the  distribution  of  signal  amplitudes.  It  is  simplest 
to  average  the  m.g.f.  in  (2.3)  and  invert  it  to  calculate  the  cumulative 
distribution  of  Z  and  thence  the  average  probability  of  detection.  Four  typical 
distr ibutions  of  the  fading  amplitudes  were  enounced  by  Swerling  [10],  and  we 
list  them  here  with  the  resulting  m.g.f.'s  of  Z: 


Case  I:  D  =  ADq,  where  A  has  a  Rayleigh  distribution, 

p( A)  =  (A/Aq2)  exp  (-A2/2Aq2),  A  >  0,  (2.6) 

yields  the  m.g.f. 


h  1  ( u )  =  (1  +  u)"(M“1}[1  +  (B1  +  Du] 

b1  =  y2  <d2>  =  <et>/n, 


(2.7) 


<•>  indicating  an  average  with  respect  to  the  distribution  of  signal  amplitudes. 


Case  II:  dj  -  Ajd0,  V j ,  where  the  Aj  are  statistically  independent  and  have 


Rayleigh  distributions  like  that  in  (2.6),  yields  the  m.g.f. 


I, 


h2(u)  =  [1  +  0  ♦  B2)u]  , 

B„  =  %  <d  ,2>  =  <E  .  >/N  =  B,  /M. 


(2.8) 


The  distribution  of  the  statistic  Z  is  a  scaled  chi-squared  distribution  with  2M 
degrees  of  freedom. 


Case  III:  D  =  ADq,  where  A  has  the  density  function 


p( A)  =  (A3/2Aq2)  exp  (-A2/2A02),  A  >  0, 


(2.9) 


yields  the  m.g.f. 


h  (u)  =  (1  +  u)'(M_2)[1  ♦  (1  ♦  B3)u]“2, 


B„  =  %  <D  >  =  <E_>/2N. 


(2.10) 


^  '  T 


Case  IV:  dj  =  Aj dg »  V j ,  where  the  Aj  are  statistically  independent  and  have  the 
distribution  of  Case  III,  yields  the  m.g.f. 


hjj(u)  =  (1  +  u)M[l  ♦  (1  +  Bjj)u]_2M, 


B„  =  y  <d.  >  =  <E . >/2N  =  B_/M. 


(2.11  ) 


The  distributions  of  the  sum  Z  arising  from  these  "Swerling  cases"  are  given 
by  di  Franco  and  Rubin  [11,  ch.  11].  For  cases  III  and  IV  the  cumulative 
distributions  are  quite  complicated,  involving  hypergeometric  functions  of  orders 
depending  on  M.  For  case  IV,  for  instance,  calculating  the  detection  probability 
Qd  requires  summing  M  terms,  each  containing  an  incomplete  gamma-function  of 
order  M+k-1,0<kSM. 

In  general,  if 


,-v 


is  the  m.g.f.  of  the  positive-valued  random  variable  S,  representing  the  total 
strength  of  the  fading  echoes,  the  m.g.f.  of  the  random  variable  Z  in  (2.1)  is, 
by  (2.3), 

h(u)  =  (1  +  u)~M  G(-p~)  •  (2.13) 

Its  singularities  will  always  lie  in  the  left  halfplane.  In  particular,  if  the 
signal  strength  has  a  gamma  distribution,  of  which  Swerling's  cases  I  -  IV  are 
special  instances, 

p(S)  =  [r(k)]”1  (k/S)k  Sk_1  e"kS/S,  (2.14) 

where  3  is  the  mean  total  signal-to-noise  ratio,  then 

G(u)  =  (1  +  frk  (2.15) 

and  the  m.g.f.  of  Z  is 

h(u)  =  (1  +  u)k_M  (1  +  bu)“k,  b  =  1  +  S/k.  (2.16) 

For  cases  I  -  IV,  k  =  1 ,  M,  2,  2M  respectively. 

The  method  of  integration  along  a  parabolic  contour,  as  in  (1.21),  has  been 
successfully  applied  to  evaluating  the  tail  probabilities  for  both  unfading 
signals  and  signals  having  random  strengths  described  by  distributions  p(S)  of 
the  form  in  (2.14).  A  bound  has  been  derived  for  the  truncation  error  incurred 
by  cutting  off  the  integration  in  (1.21)  at  a  finite  value  of  y  [12].  In  most 
instances  it  suffices  to  stop  the  numerical  integration  when  the  absolute  value 
of  the  integrand  falls  below  a  fraction  e  of  the  accumulated  sum  times  the  step 
size  Ay,  whereupon  the  relative  truncation  error  will  be  less  than  e.  The  number 


of  steps  of  numerical  integration  required  to  attain  a  certain  accuracy  is  more 


or  less  independent  of  the  number  M  of  signals  processed  by  the  radar  receiver. 

(b)  Detection  in  Noise  of  Unknown  Level 

When  the  radar  is  being  jammed  by  a  transmitter  of  Gaussian  noise  spread 
over  a  frequency  band  much  wider  than  that  of  the  target  echoes,  the  receiver  is 
faced  with  the  problem  of  detecting  them  in  noise  of  unknown  and  random  spectral 
density.  A  stratagem  that  is  often  adopted  utilizes  K  samples  r'  of  the 
quadratically  rectified  output  of  a  narrowband  pass  filter  tuned  to  the  same 
frequency  as  that  of  the  echoes,  or  to  a  nearby  frequency;  these  samples  are 
taken  at  times  when  no  signals  are  expected  to  be  present,  and  it  is  presumed 
that  they  have  the  same  probability  distribution  as  the  noise  components  of  the 
terms  r  j  in  (2.1).  The  strength  of  that  noise  is  then,  within  an  inessential 
constant  of  proportionality ,  estimated  by  the  value  of  the  statistic 
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When  these  samples  rj  indeed  represent  noise  alone,  the  statistic  Z’  has  a 
scaled  chi-squared  distribution  with  2K  degrees  of  freedom  and  is  independent  of 
Z  in  (2.1).  The  decision  level  with  which  the  statistic  Z  is  compared  in  order 
to  decide  whether  any  radar  echo  is  present  is  made  proportional  to  Z',  or 
equivalently,  the  receiver  forms  the  statistic 

X  =  Z  -  bz\  (2.18) 

and  it  decides  that  a  signal  is  present  whenever  X  >  0.  The  constant  6  is 
selected  to  achieve  a  pre-assigned  false-alarm  probability 


Qq( B)  =  Pr  (X  >  0 


(2.19) 


and  one  wishes  to  calculate  the  probability 


Qd(B;  S)  =  Pr  (X  >  0|  H  )  (2.20) 

of  detecting  a  sequence  of  M  radar  echoes  of  total  signal-to-noise  ratio  S. 

If  the  signal  to  be  detected  has  a  fixed,  known  strength,  as  in  (2.3),  the 
statistic  X  has  the  m.g.f. 

n(u)  =  E(e  UX)  =  (1  +  u)  M(1  -  Bu)  K  exp  (-  (2.21) 

with  S  the  signal-to-noise  ratio  defined  in  (2.3).  The  probability  of  detection 
is  then  given  as  in  (1.10)  by 


Qh(B  ; 
d 


S)  = 


(1 


,-M,. 

u)  (1  - 


Bu)  K  exp  (- 


Su 


_)  du 

u  2-iri  * 


(2.22) 


where  C_  is  a  contour  parallel  to  the  imaginary  u-axis  and  passing  between  the 
origin  and  the  essential  singularity  of  the  integrand  at  u  =  -1 .  One  evaluates 
(2.22)  when  the  signal  strength  S  is  so  small  that  the  expected  value  of 
X  =  Z  — B Z '  is  negative;  here 


E(X)  =  S  +  M  -  BK. 


(2.23) 


When  E(X)  >  0,  on  the  other  hand,  one  evaluates 


+ 


(1+  u)~M(1 


-  Bu)  K  exp  (- 


Su 


du 

2iri 


(2.24) 


along  a  contour  lying  in  0  <  Re  u  <  1/6. 

The  prescription,  "Choose  the  alternative  hypothesis  Hi  when  X  = 

Z  -  BZ'  >  0,"  is  equivalent  to  the  F-test  of  the  analysis  of  variance,  and  the 
F-statistic  is 
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F  =  KZ/MZ’. 


(2.25) 


The  false-alarm  probability  Qq(B)  is  related  to  the  central  F -distribution,  ana 
the  detection  probability  Qd(S;  S)  to  the  noncentral  F-distribution,  the  signal- 
to-noise  ratio  S  playing  the  role  of  the  noncentrality  parameter.  In  [13]  the 
computation  of  the  noncentral  F-distr ibution  by  integrating  (2.22)  or  (2.24) 
numerically  along  a  straight  vertical  contour  was  treated,  and  a  bound  on  the 
error  incurred  by  truncating  the  range  of  numerical  integration  was  presented. 

It  was  shown  that  the  simple  rule  that  the  summation  be  stopped  when  the  absolute 
value  of  the  integrand  falls  below  a  fraction  e  of  the  accumulated  sum  times  the 
step  size  Ay  suffices  to  bound  the  relative  error  in  the  computed  probability 
within  that  same  fraction  e. 

If  the  radar  echoes  are  fading  as  described  in  part  (a),  the  average 
detection  probability  is  obtained  by  averaging  the  integrals  in  (2.22)  and  (2.24) 
with  respect  to  S, 

Q,(B;  S)  =  -  /  u'\l  +  u)~M(1  -  Bu)”K  C(— — )  (2.26) 

Q  /  1  +  U  <:7I  1 

*X 

1  -  Qd(B;  S)  = 

where  S  is  the  average  total  signal-to-noise  ratio  and  G(u)  is  the  m.g.f.  of  the 
distribution  of  the  random  signal-to-noise  ratio  S,  as  in  (2.12).  Again  it  is 
possible  to  evaluate  these  probabilities  by  numerical  integration  along  a 
vertical  straight  line  or  other  suitable  contour.  Bounds  on  the  truncation  error 
when  the  m.g.f.  is  given  as  in  (2.15)  have  been  determined  [12]. 


u  1  (1 


u)'M(l  - 


Bu)"K  G[ 


du 
2ni  ’ 


(2.30) 


3.  Robust  Radar  Detectors  (J.  A.  Ritcey) 


Mean-level  detectors  (MLD)  are  commonly  used  in  radar  to  maintain  a  constant 
false-alarm  rate  (CFAR)  when  the  background  noise  level  is  unknown.  In  Section  2 
we  described  the  application  of  saddlepoint  integration  techniques  to  radar 
detection  problems  in  which  the  test  statistic  is  compared  with  either  a  fixed 
threshold  or  with  an  adaptive  threshold  composed  of  a  weighted  estimate  of  the 
mean  level  of  the  background  noise,  derived  from  a  given  reference  channel 
[12].  The  method  was  developed  for  an  arbitrary  number  of  pulses  noncoherently 
integrated  and  for  both  a  nonfluctuating  and  a  chi-square  fluctuating  target.  We 
have  also  considered  a  variant  of  the  basic  MLD  that  is  more  robust  in  the 
presence  of  interference  in  the  reference  channel,  and  we  have  obtained  exact 
analytical  results  for  the  probabilities  of  detection  and  false  alarm  when  the 
interference-to-noise  ratio  (INR)  is  nonzero.  In  practice,  the  interference 
could  result  either  from  returns  from  real  objects  such  as  other  targets  or 
clutter,  or  from  pulsed-noise  jamming.  Here  we  consider  only  Swerling  II  target 
returns  and  only  single-pulse  processing  (no  noncoherent  integration).  The 
approach  is  to  write  the  detection  probability  as  a  contour  integral,  which  is 
evaluated  by  the  method  of  residues.  Variants  of  this  system  in  which 
noncoherent  integration  is  employed  or  other  target  models  are  considered  would 
require  numerical  contour  integrations.  In  the  model  we  consider  here, 
developing  an  analytical  expression  for  the  moment-generating  function  (m.g.f.) 
of  the  reference-channel  estimate  is  the  primary  mathematical  problem. 

Radar  detection  of  a  known  signal  in  additive  white  Gaussian  noise  of 
unknown  variance  is  often  accomplished  by  comparing  the  test  statistic  for  a 
single  range  cell  with  an  adaptive  threshold  equal  to  a  scaled  estimate  of  the 
unknown  noise  variance.  For  a  system  that  quadratically  rectifies  the  output  of 
a  matched  filter  to  obtain  the  test  statistic,  the  problem  can  be  modeled  by  the 
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following  hypothesis  testing  problem: 

Hq:  v  =  ^  (x  +  y  } 

H1 :  v  «  j  ((x+s)2  +  (y  +  t)2);  j  d  J 2  =  s2  +  t2.  (3.1) 

The  random  variables  x,  y  are  independent  Gaussian  noise  samples  with  zero  mean 
and  unit  variance.  The  target  echo  provides  a  fixed  signal-to-noise  ratio  (SNR) 
for  a  steady  target,  and  the  amplitude  can  be  taken  as  a  random  variable  when  the 
echo  fluctuates.  The  phase  of  the  complex  amplitude  d  =  s  +  it  is  uniformly 
distributed  over  (0,  2it ) ,  and  the  total  SNR  is 

S  =  | d  | 2 /2 .  (3.2) 


Implementing  the  generalized  likelihood  ratio  test,  the  system  decides  for 
hypothesis  Hq  ("signal  absent")  or  for  hypothesis  H1  ("signal  present")  according 
to  whether 


< 


ar 


(3-3) 


where  r  is  an  estimate  of  the  noise  variance  and  a  is  a  positive  thresnoid  param¬ 
eter.  Typically  the  radar  uses  the  n  range  cells  surrounding  the  cell  under  test 
to  compute  an  estimate 


n 


j=1 


where 


and  the  x j ,  yj ,  j  =  1,2,  ...,  n,  are  independent  and  identically  distributed 
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(i.i.d.)  zero-mean  Gaussian  random  variables;  we  have  incorporated  a  scaling 


factor  n  1  into  the  threshold  parameter  a.  In  a  locally  homogeneous  noise 

environment  the  x.p  y.  all  have  common  unit  variance.  The  performance  is 
J  J 

completely  specified  by  the  false  alarm  and  detection  probabilities 


(3.6) 


(3.7) 


where  f(r)  is  the  probability  density  function  (p.d.f.)  of  the  estimate  r  and 
p^(v),  i  =  0,  1 ,  is  the  p.d.f.  of  v  under  each  hypothesis.  An  equivalent 
hypothesis  test  is 

H1 

z  =  v  -  ar  ^  0  (3.6) 

H0 


where  the  random  variable  z  has  a  p.d.f.  p  (z)  whose  Laplace  transform  is  given 
by 


h  (u) 
z 


•/ 


-uz 


Pz(z)dz  = 


/ 


-uv 


p( v  )dv 


/ 


ar  u 


f  (r  )  dr 


(3.9) 


or 


hz(u)  =  <e  UZ>  =  <e  uvXeuar>  =  h(u)  d(-au)  (3.10) 

where  h(u),  d(u)  are  the  moment-generating  functions  (m.g.f.)  of  v  and  r 
respectively.  If  we  specify  a  Swerling  II  target  model,  the  m.g.f.  is  given  by 


-I  2 

and  S  =  -  <|d[  >  is  the  average  SNR.  The  estimate  r,  in  a  homogeneous 
environment,  has  the  m.g.f.  d(u)  given  by 


d(u)  =  (1+u)  , 


(3.12) 


The  detection  probability  can  be  written  as  a  contour  integral  involving 
these  m.g.f. 's.  Since  h  (u)  and  P2(z)  are  related  by  an  inverse  Laplace 
transform,  we  find  directly  from  (3.8)  that 


or 


(3.13) 


Q 


d 


u  1  h(u)  d(-au) 

2iti 


(3.14) 


where  the  contour  of  integration  C_  consists  of  a  vertical  path  in  the  complex 
u-plane  crossing  the  negative  real  axis  at  u  =  c^ ,-b  1  <  c1  <  0.  The  contour  is 

closed  in  an  infinite  semi-circle  in  the  left  half-plane  (LHP).  Since  d(u)  is 
the  m.g.f.  of  a  positive  random  variable,  the  only  singularities  lie  in  the 
LHP.  Therefore,  d(-au)  is  analytic  in  the  LHP  and  (3-14)  can  be  immediately 
evaluated  in  terms  of  the  residue  at  the  simple  pole,  u  =  -b  1 .  Thus,  the 
detection  probability  is  simply  related  to  the  m.g.f.  of  the  estimate  r  by 


Qd  =  d(a/b);  b  =  1  +  S,  (3.15) 

and  the  false-alarm  probability  is 

Qq  =  d(a ). 

The  key  point  of  this  result  is  that  given  a  closed-form  analytical  expression 
for  the  m.g.f.  d(u)  we  immediately  have  a  result  for  the  detection  performance. 
We  are  encouraged,  therefore,  to  investigate  other  systems  that  develop  a 
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different  estimate  r,  but  one  whose  m.g.f.  can  still  be  obtained. 

To  reduce  the  degradation  caused  by  interfering  targets  falling  within  the 
reference  channel,  Rickard  and  Dillard  [14]  proposed  a  class  of  detectors  S1*, 

£  =  0,  1,  2 . n-1 ,  where  the  l  largest  range  cells  are  censored  from  the 

noise-variance  estimate.  If  the  quadratically  rectified  output  for  the  j-th 
range  cell  is  denoted  by  Zj  ,  the  detector 2>k  uses  the  unsealed  estimate 


r  z< 


(3.16) 


where  k  =  n-1  and  0  £  z,.,  £  z,,.,  £  ...  £  z,  .  £  ®  are  the  order  statistics  of 

( 1  )  (2)  (n ) 

the  sample  (z^  ,  z 2>  ...,  z  ).  Since  r  is  a  biased  estimator  of  the  unknown  noise 
variance  we  consider  the  more  general  statistic 


r  2u>  *  0ZM 


(3.17) 


for  some  constant  c  >  0.  The  contour  integral  (3.14)  requires  the  m.g.f.  of  r, 


For  a  homogeneous  noise  environment  where  the  Zj's  are  i.i.d.  with  p.d.f. 
p(z)  =  e  z,  z  £  0,  the  m.g.f.  of  r  can  be  shown  to  be  given  by 


d(u>  -  ♦  (££j)  u]‘ 


(3.18) 


From  (3.18)  we  note  that  the  noise-level  estimate  will  be  unbiased  if  we  scal< 
by  k”1  and  choose  c  =  n  +  1  -  k.  Then 


§  ’  i  Ezu>  *  ln*,‘k)z(k)l 


(3-19) 


.  •  W  «.*  »  *  _4  W* 


is  an  unbiased  estimate  of  the  unknown  noise  variance  that  is  robust  with  respect 


to  the  presence  of  several  large  interference  pulses  in  the  reference  samples 

(V  znK 

The  performance  of  the  CMLD  in  a  multiple-target  environment  can  be  obtained 
directly  from  (3.15),  given  the  m.g.f.  d(u).  It  can  be  shown  that  the  joint 

m.g  .f . 

h£(u,s)  =  <exp[-u(z^1)  +  ...  +  j  )-sz(j< )  ]>,  (3.20) 


where  z  S  z(J)  S 


<  z  < 
■  z(k)  “ 


'(n) 


are  the  order  statistics  drawn  from 


an  independent  and  exponentially  distributed  sample  (z^,  z 2>  ....  z^  ...,  z  )  and 


where 


is  given  by 


where 


<z . >  =  » 

1  I 


(. 


-1 


i  =  1 ,  ....  I, 
i  =  £+1 ,  ....  n , 


a  =  1  +  INR, 


k  k 

EL-  E  Ei 


Vu< 


i=0 


k  k 


1  + 


(a-1 ) (i+r . (i ))  +  0  .  -1 
_ J _ ij 


V1  V1  Vi  =  1  J  =  1 


£1  *  £2  *  •••  *  ££-i 


iJj  =  n  +1  -  j, 

'J 


0  =  s  +  u  (k  -  j ) , 


:3.2D 


(3.22) 


The  multiple  sum  in  the  second  line  of  (3.22)  is  present  only  when  i  <  £.  The 
reason  we  have  evaluated  the  joint  m.g.f.  h^(u,s)  is  twofold:  First,  the  CMLD 
with  arbitrary  bias  constant  c  has  a  reference  channel  m.g.f.  given  by  h^(u,cu) 
when  £  inter ferers  are  present.  Second,  we  can  evaluate  the  performance  of  the 
so-called  Order  Statistic  Detector  (OSD)  described  by  Rohling  [15],  which  was  the 
estimate  r  =  Z(k),  a  single  order  statistic.  This  system  may  be  considerably 
easier  to  implement  than  the  CMLD. 

As  an  example  of  a  particular  result.  Figure  3.1  shows  the  additional  SNR 
required  (often  called  the  CFAR  loss)  for  the  MLD,  CMLD,  and  OSD  in  a  multiple 
target  environment  to  achieve  a  given  Qd  at  fixed  Qq  over  that  required  by  the 
Neyman-Pearson  detector,  which  is  the  optimal  detector  when  the  background  noise 
level  is  known.  Notice  that  although  the  CFAR  loss  increases  with  INR  for  the 
MLD,  the  CMLD  and  OSD  have  a  bounded  loss  as  the  INR  •+  ®.  These  results  have 
been  more  fully  reported  in  a  paper  to  appear  in  the  IEEE- AES  [16]  and  in  a  Ph.D. 


in  a  multiple  target 


4.  The  Integrated  Output  of  a  Linear  Rectifier 


s 

% 


* 
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In  the  receiver  of  a  radar  searching  for  a  target  at  a  particular  range,  the 
output  of  a  narrowband  filter  matched  to  the  echo  signal  is  rectified  and  sampled 
at  an  appropriate  time  to  provide  a  datum  v.  The  data  v, ,  v„,  v,.  acquired 

1  c.  M 

during  M  successive  interpulse  intervals  are  summed,  or  "integrated,"  to  yield  a 
statistic 

V  *  v1  +  v2  +  ...  ♦  vM,  (4.T) 

and  if  this  exceeds  a  certain  decision  level  Vq,  the  radar  decides  that  a  target 
is  present. 

When  the  noise  is  Gaussian  and  the  rectifier  has  a  quadratic  characteristic, 
the  detection  probability  q+(VQ)  =  Pr  (V  >  V  )  is  given  by  the  M-th-order 
Q-function  [9,  p.  219].  If  the  rectifier  is  linear,  on  the  other  hand,  the 
distribution  of  V  is  not  known  in  analytical  form.  Marcum  [18]  used  the 
Edgeworth  series  to  approximate  the  false-alarm  and  detection  probabilities  of  a 
system  integrating  the  outputs  of  a  linear  rectifier.  This  problem  is 
conveniently  handled  by  the  method  described  in  Section  1 . 

With  appropriate  normalization  we  can  write 

vk  =  [<sx  +  xk)2  +  (sy  +  yk)2]1/2  (4.2) 

where  xk  and  yk  are  independent  Gaussian  random  variables  with  mean  zero  and  unit 
1  2  2 

variance,  and  S  =  —  (s  +  s  )  is  the  input  signal-to-noise  ratio,  assumed  the 

x  y 

same  for  all  k,  M  k  ^  H.  Then  the  probability  density  function  of  the  vk’s  is 

p(v)  =  v  exp(-  ^  v2  -  S)  IQ(/2Sv),  v  >  0,  (4.3) 

where  Iq(«)  is  the  modified  Bessel  function  [9,  p.  1 66 ] .  The  moments  of  vk  were 
given  by  Rice  [19,  eq.  (3.10-12),  p.  107]  as 
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Mn  =  E(vkn)  =  2n/2r(|  +  1)  1 F1 ( -n/2 ,  1;  -S) 


(4.4) 


in  terms  of  the  confluent  hypergeometric  function 


1  F1  (a ,  1  ;  x ) 


=  1  + 


a+r)  x 


r  =  1 


U)  (ri)2* 


(4.5) 


When  S  =  0  the  m.g.f.  of  vk  is  given  in  closed  form  by  [19] 


h^iz)  =  E[exp(-vkz)]  =  1  -  (2n)172  z  eZ  72  erfc  z, 


(4.6) 


where 


erfc  z 


is  the  error-f unction  integral.  The  m.g.f.  of  V  is  then 

h(z)  =  [hQ(z)]M. 


(4.7) 


(4.  8) 


When  n  is  large,  the  n-th  moment  is  asymptotically 

v  *  2n/2T(|  ♦  De'^72  I  (/2(n+1 )  S) 
n  2  0 

by  eq.  (13*5.1 3),  p.  508  of  Abramowitz  and  Stegun  [20].  Using  Stirling's  formula 
for  the  factorial  and  the  asymptotic  formula  for  the  modified  Bessel  function,  we 
find  for  n  >>  1 

an  =  wn/n  ~  (8irnS)  n  exp[-^-  +  (2nS)  J. 


These  coefficients  decrease  so  rapidly  with  increasing  n  that  the  power  series  in 
(1.4)  for  the  m.g.f.  of  the  converges  everywhere  in  the  complex  z-plane. 

The  first  four  coefficients  of  this  series 


29 


ho(z)  =2^V_Z) 


k=0 


are 


ao  =  1  *  ai 


(n/2)  1F1  (-  |,1 ;-S) 


1/2 


The  rest  can  be  evaluated  by  the  recurrence 

2(m+1  +  S) (m-1 )a  -  ma 


m-2 


m+2 


(4.9) 


a2  =  1  +  S,  a3  =  (tt/8)  ]  F1  (-3/2, 1  ; -S) .  (4.10) 


(4.11) 


(m  -  1 ) (m  +  2) 

so  that  the  series  in  (4.5)  needs  to  be  summed  only  twice.  One  can  then  use  the 
recurrence  in  (1.7)  to  determine  the  coefficients  of  the  series 


Sin  hQ(z)  =  ^  bL.(-z)k.  (4.12) 

k  =0 

Multiplying  these  by  M  yields  the  terms  in  the  expression  (1.17)  for  the  phase 
'Hz)  of  the  integrands  of  (1.9),  (1.10),  and  (1.12). 

The  false-alarm  probability  Qq  =  P r  (V  >  V^|  S  =  0)  can  be  computed  by 
numerical  quadrature  of  (1.21),  using  (4.6)-(4.8)  to  calculate  h(z)  on  the  path 
of  integration.  The  saddiepoint  Xq  was  calculated  by  Newton's  method  as  in 
(1.16),  and  the  curvature  <  of  the  parabolic  path  of  integration  was  determined 
by  (1.22).  By  combining  this  method  with  the  secant  method  we  were  able  to 
determine  the  decision  level  Vq  to  attain  a  pre-assigned  false-alarm  probability 
Qq.  The  detection  probability  =  Pr  (V  >  Vq J  S  >  0)  could  then  be  calculated 
by  (1  .21 ) ,  in  which  the  m.g.f .  h(z)  was  computed  as  in  (4.9)— (4.12),  with 
(4.8).  This  method  worked  very  efficiently  and  enabled  us  to  compute  the  curves 
shown  in  Fig.  4.1,  which  supplement  the  one  given  by  Marcum  [18,  p.  253 J.  These 


show  the  ratio  in  dB  of  the  signal  strength  required  when  using  a  linear 


PvT" 


< 


«  *  irk  ■ 


^  \r* 


rectifier  to  that  required  when  using  a  quadratic  rectifier  to  attain  a  given 
probability  Qd  of  detection,  for  Qd  =  0.5,  0.99,  0.999,  0.9999  and  for  various 
numbers  M  of  pulses  integrated.  The  false-alarm  number 


Nfa  =  (M  In  2)/Qq 


was  set  at  10  . 


Fig.  4 
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.  1.  Ratio  in  dB  of  SNR  required  by  linear  detector  to  that  required  by- 

quadratic  detector  to  attain  detection  probability  Qd ,  versus  number  M 
of  pulses  integrated.  Curves  are  indexed  with  the  values  of  Qd> 
False-alarm  number  N^.  =  10^. 
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5.  Distribution  of  the  Filtered  Output  of  a  Quadratic  Rectifier 
In  a  fundamental  paper  [21]  Kac  and  Siegert  showed  how  to  calculate  the 
moment-generating  function  (m.g.f.)  of  the  filtered  output  of  a  quadratic 
rectifier  whose  input  is  Gaussian  noise  and,  possibly,  a  coherent  signal.  Only 
in  one  very  special  case,  however,  did  they  invert  the  m.g.f.  to  obtain  the 
probability  distribution  of  the  output.  In  general ,  numerical  methods  are 
necessary,  and  saddlepoint  integration  has  proved  to  be  a  most  efficient  one. 

Let  the  input  to  the  rectifier  be  a  narrowband  Gaussian  random  process 

v ( t )  =  Re  [ SQ(t )  ♦  N(t)]elfit  (5.1) 

with  carrier  frequency  ft.  Here  SqU)  is  the  complex’ envelope  of  a  coherent 
signal,  and  N(t)  is  a  stationary  circular  Gaussian  random  process  with  mean 
zero  and  complex  autocovariance  function 

$(t1-t2)  =  1/2E[N(t1)N«(t2)],  (5.2) 

with  E[ N(t1 ) N ( 1 2 ) 3  =  0.  The  output 

y2|s0(t)  +  N(t )  | 2 

of  the  quadratic  rectifier  is  applied  to  a  filter  with  positive  impulse  response 
K(t)  >  0.  To  be  calculated  is  the  distribution  of  the  output  w  of  that  filter  at 
time  T, 

w  =  \ J  K(t)|Sq(T-t)  +  N(T-t)|2  dt.  (5.3) 

0 


The  m.g.f.  of  the  random  variable  w  is  [21] 


j  -  i-A  c  ;  -  i,  u  \  Ct  j  j  CA^  I  2  1  +  ^  z  J  t 

k  k 


where  the  are  the  eigenvalues  of  the  integral  equation 


( 


4>(u-v)  K (v )  gR(v)  dv  =  Akgk(u),  0  <  u  <  T, 


(5.5) 


whose  eigenfunctions  are  g>K(t),  normalized  so  that 

T 


K(t)  g  *(t)  g  ( t )  dt  =  6  . 

k  m  km 


(5.6) 


In  (5.4) 


D  (z ) 


n  °  *  v> 


(5.7) 


is  the  Fredholm  determinant  and 


s 


k 


*(t)  K (t )  SQ«(T-t)  dt. 


(5.8) 


If  the  input  is  instead  a  lowpass  Gaussian  process  sQ(t)  +  n0(t),  where 
nQ(t)  is  Gaussian  noise  with  mean  zero  and  autocovariance  function  <(>  (t  ^ -tp)  ,  the 

m.g.f.  is 

2 

h1  (z )  =  [  D(z )  ]  /z  exp  [-  y~  »  (5.9) 

k  k 

where  the  sk  are  defined  as  in  (5.8)  with  S0*(T-t)  replaced  by  s0(T-t),  and  D(z) 
is  again  the  Fredholm  determinant  (5.7)  associated  with  (5.5). 

The  m.g.f.  h(z)  in  (5.4)  can  be  inverted  by  the  residue  theorem  to  give  the 


probability  density  function  p(w)  of  the  output  w  when  S(t)  =  0,  but  when  the 


product  WT^nt  of  the  bandwidth  W  of  the  input  noise  and  the  integration  time  TinL 
of  the  output  filter  is  large,  the  resulting  series  is  difficult  to  sum 
accurately  because  the  leading  terms  are  large  and  of  alternating  sign.  When 
S(t)  *  0,  one  might  convolve  the  Rayleigh-Ri ce  distributions  that  are  the  inverse 
transforms  of  the  factors  of  (5.4),  but  it  is  difficult  to  do  so  accurately, 
particularly  in  the  tails  of  the  p.d.f.  Numerical  Fourier  transf ormation  of 
h(-ito)  is  also  inaccurate  in  the  tails  of  the  distribution. 

Because  of  the  square  root  in  (5.9),  the  residue  theorem  cannot  be  applied 

to  it.  In  determining  the  p.d.f.  of  the  average  power  of  a  mean-zero  random 

process,  a  special  case  of  (5.9),  Slepian  [22]  integrated  the  Laplace  inversion 

integral  with  (5.9)  numerically  around  cuts  on  the  negative  Re  z-axis  between  the 

-V 

branch  points  -1/A^  of  [D(z)J2.  The  method  of  saddlepoint  integration  avoids 
difficulties  such  as  these  by  taking  the  contour  of  integration  in  (1.9)  or 
(1.10)  along  a  contour  close  to  the  path  of  steepest  descent  of  the  integrand  and 
far  from  its  singularities  on  the  negative  Re  z-axis. 

When  Sn(t)  =  0  and  the  input  noise  has  a  Lorentz  spectral  density, 

*(w)  =  2y(u,2  +  u2)"1,  (5.10) 


and  autocovariance  function 

♦  (t)  =  e'ulTl  ,  (5.11) 

and  when  the  output  filter  has  an  "RC"  form 

K(t)  =  6  e"BT  U(t) ,  (5. 12) 

the  m.g.f.  of  the  filtered  output  w  of  the  quadratic  rectifier,  as  shown  by  Kac 
and  Siegert  [ 21  ] ,  is 
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-  Yiv)  I  Tn)’  -  '  ’  V  ~  L0.1J5  > 

V-1  v 

where  Iv_-|(*)  is  the  modified  Bessel  function.  By  using  a  continued  fraction  for 
the  logarithmic  derivative  of  the  Bessel  function,  the  saddlepoints  of  the 
integrands  of  (1.9)  and  (1.10)  can  be  quickly  determined,  and  the  cumulative 
distribution  q  (w)  of  w  can  be  computed  by  integrating  (1.21)  by  the  trapezoidal 
rule.  In  this  way  we  produced  the  curves  for  q~(w)  shown  in  Fig.  5.1  [23]. 

Because  good  accuracy  can  be  attained  in  these  numerical  saddlepoint 
integrations  with  only  about  a  dozen  steps,  it  becomes  feasible  to  treat  problems 
in  which  the  m.g.f.  is  not  —  as  in  (5.13)  --  known  in  closed  form,  but  must 
itself  be  computed  numerically.  To  this  end  we  have  shown  that  the  m.g.f.  h(z) 
in  ( 5 . ^ )  can  be  written  as 


h(z)  =  LD(z)]  '  exp 


K( t )[o*(t ;z*)-S*( t ) ][o(t ;z)-S(t ) ]dt } 


(5.1H) 


with  S(t)  =  SQ*(T-t)  and  with  the  Fredholm  determinant  given  by  [ 2 h 3 


/ 


D(z)  =  exp  [-  /  K(t)  r(t,  t;  z)  d t ] . 
0 


(5.15) 


Here  [23,  25] 


o(t;  z) 


z)  K(u)  S' u)  du 


(5.16) 


might  be  called  the  "causal  quasi-estimate"  of  the  signal  S(t)  =  S0*(T-t).  The 


function  r(t,  u;  z),  which  is  the  solution  of  the  integral  equation 


Fig.  5.1.  Cumulative  distribution  q  (w)  of  RC- filtered  output  of  a  quaoret . : 

rectifier  whose  input  is  narrowband  Gaussian  noise  with  mean  zero  and 
a  Lorentz  spectral  density.  Curves  are  indexed  with  the  value  of  v  = 
2u/g. 
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(5.  17) 


z)  K ( v )  <f>(v-u)  dv  =  z$(t-u), 


0  <  u  <  t  <  T, 

is,  for  K(t)  =  1  and  z  real  and  positive,  the  kernel  of  the  minimum-mean-square- 

y 

error  estimator  of  the  random  process  z2N(t)  when  it  is  observed  in  white 
Gaussian  noise  of  unit  spectral  density. 

It  is  this  connection  with  estimation  theory  that  shows  a  way  of  determining 
the  m.g.f.  h(z)  numerically  without  having  to  find  the  eigenvalues  and  the 
eigenfunctions  g(<(t)  of  the  integral  equation  (5.5).  When  the  spectral  density 
*(u)  of  the  noise  is  a  rational  function  of  the  frequency,  one  can  set  up  a 
state-space  model  of  a  linear,  time- invariant  system  that  is  driven  by  circular- 
complex  white  noise  W(t)  and  generates  the  process  N(t).  The  state  equations 
have  the  usual  form 

dX 

—  =  FX  +  GW(t ) ,  (5.  16/ 

dt 

where  X(t)  is  a  column  vector  of  n  state  variables,  with  2n  the  degree  of  the 
denominator  of  the  spectral  density  $(uj)  ,  G  is  a  constant  n-element  vector,  and  F 
is  the  nx.n  dynamical  matrix  of  the  system.  For  an  appropriate  n-element  row 
vector  Z,  the  input  noise  can  be  represented  as 

N(t)  =  CX(t).  (5.19) 

Then  it  can  be  shown  that  the  kernel  r(t,  u;  z)  is  given  by 

r(t,  u;  z)  =  Ck^ (t ,  u) ,  (5.20) 

where  the  n-element  column  vector  k^(t,  u)  obeys  the  linear  differential  equation 
9k  (t,  u) 

— -~rr -  =  [F  -  I  (t,  t  )K  (t  )C]k  (t,  u),  0  <  u  <  t  <  T,  (5.21) 


(5.22) 


lim  k  (t,  u)  =  k  (u,  u)  =  £(u;  z)C 
-z  -  z  ^  ~ 

t+u 

with  the  nxn  matrix  £(t;  z)  the  solution  of  the  matrix  Ricatti  equation 


9g(t;  z)  +  +  + 

- -  =  Fg  +  £F  -  C  K(t)C|  +  zQGG, 


|(0;  z)  =  zAq, 


(5.23) 


where  Q  is  the  spectral  density  of  the  driving  white  noise  W(t)  in  (5.18)  and 
is  the  steady-state  covariance  matrix 


Aq  =  E[X(t)  X+(t)3 

and  obeys  the  equation 

FAq  +  AqF  +  QGG+  =0.  (5.24) 

The  causal  quasi-estimate  o(t;  z)  in  (5.16)  is 


o(t;z)=Co(t;z),  (5. 25) 

where  the  n-element  column  vector  o ( t ;  z)  is  the  solution  of  the  vector 
differential  equation 


3o(t;  z) 

at 


Fo(t;  z)  +  £(t;  z  )C+K  (t )  [  S(t )  -  Co(t;  z)], 
2 ( 0 ;  z)  =  0. 


(5.26) 


With  J(T;  z)  denoting  the  integral  in  the  exponent  of  the  m.g.f.  h(z)  in  (5.14), 
and  with  E(T;  z)  =  In  D(z),  one  solves  the  differential  equations  (5.23)  and 
(5.26)  along  with 


9 E ( t ;  z) 
at 


K(t)Cg(t;  z)C  + , 


ECO;  z)  =  0, 


(5.27) 


(5.28) 


*2 


3 J( t ;  z) 


=  K ( t ) [ o  * ( t ;  z  * )  -  S#(t ) ][ o(t ;  z)  -  S(t)], 
J(0;  z)  =  0, 


over  the  range  0  <  t  <  T.  This  can  be  accomplished  by  standard  routines  for 
solving  sets  of  first-order  differential  equations.  At  t  =  T  one  has  the 
functions  E(T;  z)  and  J(T;  z)  for  use  in  the  phase 


V(z)  =  -E  ( T ;  z)  -  |  J(T;  z)  +  wz  -  In  (±z) 


(5.29) 


in  evaluating  the  saddlepoint  integral  0.21).  One  carries  out  this  procedure  at 
each  point  z  on  the  contour  of  integration  after  (1.21)  has  been  replaced  by  a 
suitable  quadrature  formula  such  as  the  trapezoidal  rule. 

The  derivatives  of  the  phase  y(z)  for  real  values  of  z  needed  for  finding 
the  saddlepoint  points  xQ+  and  xQ_  by  Newton's  method  as  in  (1.16)  and  for 
determining  the  curvature  <  of  the  path  of  integration  by  (1.22)  can  be  computed 
by  solving  along  with  (5.23),  (5.26),  (5.27),  and  (5.28)  additional  first-order 
differential  equations  obtained  by  differentiating  these  equations  and  their 
initial  conditions  with  respect  to  z.  This  program  has  been  carried  through  for 
RC-filtering  of  the  output  of  a  quadratic  rectifier  whose  input  is  narrowband 
Gaussian  noise  with  a  Lorentz  spectral  density  plus  a  signal  with  constant 
complex  envelope.  Good  agreement  was  observed  with  the  cumulative  distribution 
computed  by  using  (5.4)  for  the  m.g.f.  with  the  factor  [D(z)]_1  given  by  (5.13) 
and  with  the  constants  sk  calculated  by  (5.8)  [23]. 


6.  Distribution  of  the  Average  Power  of  a  Gaussian  Process 


A  special  case  of  the  problem  treated  in  Section  5  is  that  of  finding  the 
distribution  of  the  average  power 


/  lN(t)i2  dt 

0 


(6.1) 


of  a  narrowband  Gaussian  random  process  Re  N(t)elfit  with  carrier  frequency  and 
mean  zero,  or  of  a  lowpass  Gaussian  process  n(t), 


w' 


„-1 


T 

Ln(t)]2 


dt 


(6.2) 


again  with  mean  zero.  The  m.g.f.  of  w  in  (6.1)  is 

h(z)  =  E(e“ZW)  =  [D ( z )  ] 


with  D(z)  the  Fredholm  determinant  given  in  (5.7)  and  (5.15),  where  the  function 
r(t,  uj  z)  is  by  (5.17)  now  the  solution  of  the  integral  equation 


•*/ 


r(t,  u;  z  'J  +  r  J  r’(t>  v ;  z)  <f>(v-u)  dv  =  z<J>(t-u), 

0 

0  <  u  <  t  <  T.  (6. 3) 


For  the  random  variable  w'  in  (6.2),  the  m.g.f.  is 

h  1  ( z )  =  E(e~ZW')  =  [D(2z)]"‘/2.  (6.1) 

Slepian  [22]  calculated  the  p.d.f.  and  the  cumulative  distribution  of  w’  for 
inputs  n(t)  having  Lorentz  ("RC")  and  "RLC"  spectral  densities  by  determining 
h^(z)  analytically  and  integrating  the  Laplace  inversion  formula  numerically 
around  cuts  on  the  negative  Re  z-axis  between  pairs  of  branch  points 


zk  -  -1/Xk.  His  method  requires  finding  a  large  number  of  eigenvalues  Xk  of  the 
integral  equation  (5.5)  with  K(v)  =  2/T;  the  larger  the  time- bandwidth  product 
WT,  with  W  the  bandwidth  of  the  input  noise,  the  more  of  those  cuts  along  which 
one  needs  to  integrate.  Saddlepoint  integration  provides  a  simple  and  efficient 
way  of  inverting  the  m.g.f.  to  determine  the  distribution  of  the  average  power. 

When  the  circular  complex  Gaussian  process  N(t)  has  a  rational  spectral 
density  4>(oj)  that  can  be  put  into  the  form 


-  n,  f  U  '  +  g  U-u") 

nrm  _ m 

^ j  (w-ip  *) (w+ip  ) 

—r  m  m 


f  p  '  +  g  (w-p  ") 
m  m  m  m 

/  ,  ,  „  ,2  '2 

“t  (u>-p  " )  +  p 

m= I  m  m 


(6.5) 


where 


-ipk  =  -ip^'  +  Pk",  ipk*  =  ipk'  +  »  1  <  k  <  n, 

are  the  2n  poles  of  <t>(aj),  all  of  which  we  assume  to  be  distinct,  then  the  complex 
autocovariance  function  of  N(t)  is 


■J2  (W  exp  t>0 


T  <  0 


(6.6) 


When  this  is  substituted  into  the  integral  equation  (6.3),  the  kernel  r(t,  u;  z) 
has  the  form 


’(t,  u;  z)  c  (t)  exp  [6  (t-u)]. 


(6.7) 


in  which  the  coefficient?  6m  are  the  2n  roots  of  the  equation? 

1  +  zT  '  4>(-ip)  =  0,  p  =  B^  ,  B2>  ....  Bm>  (6.8) 

and  the  cm(t)  sati?fy  a  certain  ?et  of  linear  simultaneous  equations.  It  is  then 
possible  to  show  from  (5.15)  that  the  Fredholm  determinant  D(z)  is  given  by  [23] 

det(D2'1M*'1  )  det  -  B^C^1  §1  D,  _1  ) 

D(z)  =  - : - - - ,  (6.9) 

det  (C1  -  B2C2‘ 

where  the  n  *  n  matrices  involved  are  defined  by 


(~1  )mk  "  [(Bk  - 


(S2>.k  =  C(Bk+n  +  V)T]  • 


(Vmk  =  [(Bk  +  C)T]  * 


(SJm,  =  C ( B .  -  u  )T] 

-2  mk  k+n  m 


--1 


(D. )  ,  =  exp  ( B  T ) 6  , 

-1  mK  rri  ms' 


m  mk ' 


’  exp  <  W>  S„k’ 


(M*  )  ,  =  exp  (-u  *T)  6  ,  . 

mk  m  ml* 


(6.10) 


m  '  mk 


The  root?  Bm  of  (6.8)  are  taken  to  have  been  arranged  so  that  when  the  point 
z  moves  off  to  infinity  in  0  <  arg  z  <  -n ,  the  first  n  root?  trace  paths  finishing 
in  the  right  half-plane  and  the  rest  (k  =  n+1 ,  ....  2n)  trace  paths  finishing  in 
the  left  half-plane.  When  z  is  real,  8  =  -B  *,  1  <  m  <  n.  The  phase  H'(z)  of 

the  integrand  of  (1.21)  is  now 


H'(z)  =  -In  D(z)  +  wz  -  In  (±z). 


(6.11) 


The  form  (6.9)  is  designed  to  avoid  overflow  in  the  computer.  Because  most 
computer  libraries  contain  routines  for  multiplying  complex  matrices  and  forming 
their  determinants  and  inverses,  it  is  unnecessary  to  work  out  the  Fredholm 
determinant  in  analytic  form;  instead  one  simply  evaluates  (6.9)  at  each  point  on 


the  contour  of  integration.  As  good  accuracy  can  be  attained  in  integrating 


(1.21)  with  only  a  few  steps  of  the  trapezoidal  rule,  this  evaluation  does  not 


need  to  be  carried  out  very  often. 

The  spectral  density  of  a  lowpass  process  n(t)  has  the  simpler  form 


2u 


f 

m  m 


(6. 12) 


in  which  the  um  are  real  or  occur  in  complex-conjugate  pairs.  In  calculating  the 
distribution  of  the  average  power  w*  of  such  a  process,  as  defined  in  (6.2),  the 
form  in  (6.9)  simplifies  to 


det  (DM_1)  det  (C  -  BD~1c“1BD  1 ) 


D(z)  = 


det  (C  -  BC~ 1 B) 


(6.13) 


where 


C.  =  C(B.  -  u  )  T]  \  B  =  [(B  +  u  )  T]_ 1 , 

mk  k  m  mk  k  m 


Dmk  ■  exp  <l,kT)  W  "mk  '  “-P  (l,kT)  6mk ' 


(6.14) 


and  the  Bk  are  the  n  roots  of 


1  +  2zT  1  $>(-  ip)  =  0 


(6.  15) 


having  non-negative  real  parts.  The  phase  of  the  integrand  in  (1.21)  is  now 

$(z)  =  In  D(z)  -  w’z  -  In  (±z).  (6.16) 

In  [23]  this  method  was  applied  to  calculating  the  cumulative  distribution  of  the 
average  power  of  lowpass  Gaussian  processes  having  spectral  densities  char¬ 
acterizing  RLC  noise  as  well  as  those  wih  second-  and  fourth-order  Chebyshev 
approximations  to  a  bandlimited  spectral  density.  The  distribution  of  these  last 

two  were  compared  with  that  arising  from  a  strictly  bandlimited  spectral  density, 

-V 

as  calculated  by  (1.21)  with  h(z)  given  by  [D(z)j  2  and  D(z)  by  15.7),  in  which 
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the  eigenvalues  are  proportional  to  those  tabulated  by  Slepiao  and  Sonnets!  ;o 
for  the  prolate-spheroidal  wave  functions  [26J.  The  distribution  for  the  fourth 
order  Chebyshev  approximation,  calculated  witn  the  help  of  (6. 13)  and  (1.21), 
agreed  very  closely  with  this  one. 


7.  Intersymbol  and  Co-Channel  Interference 
In  a  typical  binary  channel  suffering  from  intersymbol  interference  the 
output  of  the  receiver  has  the  form 

y (t )  =  s (t)  +  n (t )  (7.1) 

at  time  t,  where  the  signal  s(t)  is 


s(t ) 


bka(t"kT) 

k  =  -°> 


(7.2) 


and  n(t)  is  Gaussian  noise  with  mean  zero.  The  bk's  are  +1  or  -1  with  equal 
probability  and  express  the  transmitted  binary  message;  we  assume  them 
uncorrelated.  The  output  of  the  receiver  from  a  single  transmitted  pulse  is 
a(t),  and  T  is  the  interval  between  pulses.  The  output  is  sampled  at  a  time  t  to 
produce  a  variate 


on  the  basis  of  which  the  receiver  decides  whether  a  particular  one  of  the  bk's, 
say  bg,  is  +  1  or  -1 ,  choosing  +1  if  y  >  0  and  -1  if  y  <  0.  The  probability  Pp 
of  error  is  then  the  probability  that  y  <  0  when  bg  =  +1 , 


Pg  =  Pr  (y  =  n  +  aQ  +  n  <  0 |  bQ  =  +1),  (7.4) 

p 

where  n  is  a  Gaussian  variate  with  mean  zero  and  variance  o  resulting  from  the 
noise,  and 
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tne  prime  indicating  omission  of  the  term  k  =  0. 


Calculating  the  error  probability  Pg  exactly  is  most  difficult,  for  all 
possible  combinations  of  the  "digits"  bk  =  ±1  must  be  taken  into  account, 

1  <  jkj  <  ®.  Even  if  one  assumes,  as  we  for  the  most  part  have  assumed,  that 
only  a  finite  number  n  of  pulses  before  and  after  aQ  need  to  be  included,  exact 
computation  of  Pg  involves  accounting  for  22N  combinations,  and  for  accuracy  N 
must  ordinarily  be  very  large.  Some  dozen  papers  have  appeared  in  whch  bounds  on 
the  error  probability  have  been  derived.  The  closer  the  bounds,  the  more 
complicated  they  are  to  evaluate. 

The  error  probability  can  be  calculated  quite  simply  by  saddlepoint 
integration  of  (1.9),  in  which  now  V  =  0  and  the  m.g.f.  h(z)  is  given  by 

N  , 

h(z)  =  I  I  cosh(akz)  exp(*/2o2z2  -  aQz)  ,  (7.6) 

k  =-N 

the  prime  indicating  omission  of  the  term  k  *  0  [ " 7 J .  Because  of  the  Gaussian 

p  p 

factor  exp  l/2o  s  ,  the  magnitude  of  the  integrand  of  (1.9)  drops  off  to  zero 
rapidly  along  a  vertical  contour  Re  z  =  c  >  0  when  that  contour  passes  through 
the  saddlepoint  c  =  xQ  of  the  integrand. 

In  this  problem  the  m.g.f.  h(z)  and  hence  the  integrand  of  (1.9)  possesses 

an  infinite  number  of  zeros  on  the  Im  z-axis,  lying  at  points  z  =  iy  such  that 

?  2 

aky  =  2m t  for  some  integers  n  and  k.  When  the  signal -to-noise  ratio  Sq  /o  is 
low,  the  path  of  steepest  descent  from  the  saddlepoint  Xq  on  the  Re  z-axis  runs 
into  one  of  those  zeros  and  then  passes  from  one  zero  to  another  on  the  way  out 
to  +i®.  For  large  signal -to-noi se  ratio  the  path  is  slightly  curved,  but  runs 
more  or  less  parallel  to  the  Im  z-axis.  Neither  of  these  paths  would  be  simple 
to  approximate,  and  we  have  therefore  adopted  the  original  straight  vertical  path 
z  =  xn  +  iy  through  the  saddlepoint  and  found  it  adequate. 


rt&cM fiOc? 


[i 


The  saddlepoint  Xq  is  the  solution  of  the  equation 


f’(z)  = 


_ i 

\  '  2  -1 
^  tanhCa^z)  +  0  z  -  aQ  -  z  =0, 


(7.7) 


k  =  -« 


i 


< 

K\ 


K. 


I 


H 


which  is  easily  solved  by  Newton's  method  (1.16).  As  the  saddlepoint  Xq  is  not 
needed  to  great  accuracy,  N  in  (7.7)  can  be  replaced  by  a  smaller  value  N"  such 
that  for  |k]  >  N"  the  values  of  the  signal  samples  a^  have  fallen  below  a  few 
percent  of  aQ.  The  integral  in  (1.9)  is  evaluated  by  the  trapezoidal  rule,  and  a 
simple  bound  on  the  error  incurred  has  been  determined  [27].  The  length  of  the 
computation  of  the  error  probability  by  this  method  is  now  roughly  linear  in  N, 
in  contrast  to  its  exponential  dependence  on  N  when  the  combinatorial  approach  is 
taken. 

If  all  possible  interfering  symbols  were  taken  into  account,  the  phase  of 
the  integrand  in  (1.9)  would  be 

CO 

I 

'Kz)  ^  ^  In  cosh  (a^z)  +  V2o2z2  -  aQz  -  In  z,  (7.8) 

k=-°° 

and  this  can  be  approximated  by 

N  ’ 

f ( z )  =  In  cosh  (a^z)  +  z2  -  aQz  -  In  z,  (7.9) 

k  =  -N 

where 


-N-1  » 

»,2  ■  ->o2  *  11,  \2  *  X)  ak2-  "• 

k  =  -«  k=N  +  1 

When  the  sums  in  (7.10)  can  be  evaluated,  perhaps  approximately,  the  pulses 
neglected  in  (7.6)  can  be  roughly  taken  into  account  as  contributing  to  the 
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Gaussian  noise. 

A  similar  problem  is  to  calculate  the  effect  on  a  binary  communication 
system  of  signals  leaking  from  adjacent  channels,  so-called  "cochannel 
interference."  Now  (7.5)  is  replaced  by 

L 

n  =  'y  ]  rk  cos  6k  (7.11) 

k  =  l 

where  the  rk  are  specified  positive  amplitudes  and  the  angles  0k  are 
independently  random  and  uniformly  distributed  over  (0,  2it);  L  is  the  number  of 
interfering  channels.  The  m.g.f.  h(z)  in  (1.9)  is  now 

L 

h (z )  =  |  |  IQ(rkz)  exp(‘4o2z2  -  aQz), 
k=1 

where  I0(*)  is  the  modified  Bessel  function.  The  resulting  error  probability  Pp 
can  again  be  calculated  by  integrating  (1.9)  along  a  straight  vertical  contour 
passing  through  the  saddlepoi.nt  c  =  Xq  >  0.  A  bound  on  the  truncation  error  fo.” 
this  integration  has  also  been  determined  [27]. 


8.  The  Distribution  of  Shot  and  Clutter  Noise 


Radar  clutter  noise  can  be  considered  a  random  process  of  the  form 


v(t>  Vlt_V 

k  =-® 


(8.1  ) 


in  which  the  terms  akf(t-Tk)  represent  echoes  of  the  transmitted  radar  pulses 
f(t)  from  scatterers  with  various  cross-sections  located  at  randomly  distributed 
ranges.  The  arrival  times  Tk  of  the  echoes  can  be  taken  as  a  Poisson  point 
process,  to  which  we  assign  a  constant  rate  v.  Only  for  a  finite  range  of 
indices  k  are  the  amplitudes  ak  significantly  large,  for  the  transmitting  radar 
antenna  illuminates  only  a  limited  area  of  the  scattering  surface.  We  can 
assume,  however,  that  the  interval  of  time  during  which  the  echoes  arrive  is  so 
much  longer  than  the  duration  of  each  pulse  that  the  summation  in  (8.1)  can  be 
taken  from  -®>  to  The  amplitudes  ak  in  clutter  noise  can  be  assumed  to  be 
independent  random  variables  with  a  common  probability  density  function.  We 
first  consider  the  special  case  of  shot  noise. 

(a)  Shot  Noise 

When  the  amplitudes  ak  are  all  equal  to  1,  (8.1)  represents  shot  noise. 

Rice  [19,  Sec.  1.4]  showed  that  the  m.g.f.  of  the  random  variable  v  =  v(t)  at  any 
time  t  is  then 


h(z)  =  exp 


[exp  [zf(t)]  -  1}  dt } 


(8.2) 


He  also  showed  how  to  approximate  the  distribution  of  the  amplitude  v  by 
Edgeworth's  series.  Gilbert  and  Poliak  [28]  derived  from  (8.2)  an  integral 
equation  for  the  cumulative  distribution  of  v  and  solved  it  for  some  simple  pulse 
shapes.  Lugannani  and  Rice  [29]  described  how  to  evaluate  the  cumulative  distri- 


bution  by  numerically  integrating  (1.9)  or  (1.10)  along  the  Im  z-axis  itself,  and 
they  also  gave  saddlepoint  approximations  and  Chernoff  bounds  on  this  distri¬ 
bution.  Yue  et  al .  [ 30 ]  developed  series  approximations  for  the  distribution 
valid  especially  for  low  values  of  vT,  where  T  is  the  pulse  length. 

In  a  preliminary  investigation  of  the  applicability  of  numerical  saddlepoint 
integration  to  computing  the  cumulative  distribution  of  the  amplitude  v  of  shot 
noise,  we  considered  a  pulse  containing  an  arbitrary  number  M  of  cycles  of  a 
sinewave , 

f(t)  =  cos  £lt,  0  <  t  <  T  =  2-irM/fi, 

=0,  t  <  0,  t  >  T,  (8.3) 

for  which  the  m.g.f.  from  (8.2)  is 

h(z)  =  exp  {vT[Iq(z)  -  1]},  (8.9) 

where  Iq(*)  is  the  modified  Bessel  function.  As  the  p.d.f.  of  v  is  an  even 
function,  we  need  to  consider  only  negative  values  of  the  level  V,  for  which  we 
tawe  the  saddlepoint  xQ  or.  the  positive  real  axis. 

The  m.g.f.  's  of  shot  noise,  by  (8.2),  are  entire  functions  and  have  no 
sir.g  ila-it  »es  ar/w r.ere  in  the  finite  portion  of  the  z-plane.  At  infinity, 
n-.wever,  tney  possess  numerous  singularities,  as  one  can  see  by  considering  the 

isvr,,  •.  '  :  c  of  :  )  for  j  z  j  large,  but  with  z  going  off  to  infinity  in  a 

d  i-e  :  1 1  ->n  ctner  tr.an  that  of  the  Im  z-axis.  The  phase  t(z)  of  the  integrand  in 

.s  r.ow,  for  \z\  >>  1,  Re  z  >  0,  V  <  0, 

_!/  7 

H\z)  -  (2uz)  2  vT  e  -  j  V | z  -  In  z,  (8.5) 

and  if  we  put  z  -  rei9,  the  imaginary  part  becomes 

Im  ¥(z)  -  vT(2irr)  ^  ex  sin(y  -  |-)  -  j  V  |  y  -  6.  (8.6) 
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This  must  remain  constant  on  the  path  of  steepest  descent  as  z  runs  out  to 
infinity  in  the  right  halfplane.  The  term  with  the  factor  ex  will  grow 
exponentially  unless  sin(y  -  6/2)  goes  to  zero.  This  means  that  z  must  go  off  to 
infinity  along  a  line  on  which  y  equals  a  multiple  of  tt;  along  such  a  line  the 
angle  6  goes  to  zero  as  Re  z  +  ®.  The  corresponding  term  in  Re  v(z)  is 
proportional  to 

ex  cos  (y  -  6/2), 

and  if  the  integrand  is  to  go  to  zero,  y  must  approach  an  odd  multiple  of  it.  The 
path  of  steepest  descent  thus  consists  of  hairpin  curves  opening  toward  +®  and 
approaching  the  lines  Im  z  =  ( 2k-*- 1 ) -it  asymptotically  for  all  integers  k.  On  the 
main  branch  lying  symmetrically  about  the  Re  z-axis,  k  =  0  and  1m  y(z)  =  0.  Each 
of  these  branches  possesses  one  saddlepoint;  for  k  =  0  it  is  the  point  we  have 
denoted  by  xQ. 

We  have  integrated  (1.9)  with  the  m.g.f.  in  ( 8 .  h )  along  the  main  branch  of 
the  path  of  steepest  descent  and  along  a  number  nof,^  of  the  branches  above  the 
Re  z-axis,  and  the  results  are  exhibited  in  Table  8.1.  The  branches  lying  below 
the  Re  z-axis  are  mirror  images  of  those  above  it  and  contribute  equally  to  the 
cumulative  probability  q-(V).  The  branches  were  traced  by  the  method  described 
in  Sec.  I,  part  (b),  and  the  integration  along  each  was  carried  out  by  means  of 
the  five-point  quadrature  formula  (1.27)  of  Birkhoff  and  Young  [7].  The 
contribution  of  each  branch  was  also  estimated  by  the  saddlepoint  approximation 
(1.29),  which  agreed  within  a  few  percent  with  the  result  of  the  numerical 
integration. 

The  column  marked  "Ratio"  in  Table  8.1  lists  the  ratio  of  the  contribution 
of  the  off-axis  branches  to  the  total  probability  q”(V).  For  vT  =  1 0  the  off- 
axis  branches  contributed  only  a  small  fraction  of  the  probability,  and  five  of 


Table  8.1 


q~(V):  Shot  noise 

,  sinusoidal 

pulses 

Main  branch 

Off-axis 

Ratio 

Total 

branches 

probability 

q“(  V) 

vT  =  10, 

noff  =  5 

0.32431 

-5.0l63(-5) 

-1 . 547 (-4) 

0.32426 

0 . 5826  3 ( ~ 1 ) 

4.71 56 ( —6 ) 

8.093C-5) 

0 .58268 (-1 ) 

0.31845C-3) 

9 . 4 1 86 ( —9 ) 

2.958C-5) 

0.31846 (-3) 

vT  =  1, 

nof f  =  12 

0.40362 

-0.11080 

-0.3784 

0.29283 

0.60583 ( -1 ) 

-0. 14782(-1 ) 

-0.3227 

0.4580K-1) 

0.35009 (-2) 

-0.17  276 ( -3) 

-0.05191 

0.33281 (-2) 

LRY 


0.3242 
0  .  5  8  27  (-1 ) 
0.3185C— 3) 

0.2920 
0.  -4607  ( -1  ) 
0.3320 (-2) 


them  sufficed.  For  vT  =  1 ,  on  the  other  hand,  the  off-axis  branches  contributed 


a  substantial  fraction  of  the  total  probability,  and  these  amounts  decreased  very 
slowly  and  in  an  oscillatory  fashion  as  one  moved  to  branches  farther  and  farther 
from  the  Re  z-axis.  Twelve  branches  above  the  axis  were  taken  into  account  for 
vT  =  1 .  The  column  marked  "LRY"  lists  the  values  of  q~(V)  tabulated  by 
Lugannani,  Rice,  and  Yue  [31].  The  agreement  is  excellent  at  vT  =  10,  but  not  so 
good  at  vT  =  1 .  The  smaller  the  value  of  vT,  the  more  slowly  the  m.g.f.  h(z) 
decreases  as  the  value  of  z  moves  toward  +i®. 

(b)  Clutter  Noise 

When  as  for  clutter  noise  the  amplitudes  ak  in  (8.1)  are  independent  and 
possess  a  common  probability  density  function  p(a),  the  m.g.f.  of  the  random 
variable  v  =  v(t)  is  [32] 


h(z) 


■/ 


=  exp  |v  7  [H(zf (t) )  -  1]  dt}, 


(8.7) 


where 


—  CO 


(8.8) 


is  the  m.g.f.  of  the  amplitudes  ak .  If,  as  we  shall  assume,  these  have  identical 

2 

Gaussian  distributions  with  mean  zero  and  variance  o  , 

H ( z )  =  exp  (l/2o2z2)  ,  (8.9) 


then 


-  exp  |  v 


f 


[exp(‘4o2z2[f  ( t ) ] 2 )  -  1  ] dt } . 


:  8 .  i  o ) 


We  considered  two  types  of  signal:  (i)  a  rectangular  pulse, 

r(t)  =1,  0  <  t  <  T;  f(t)  =0,  t  <  0,  t  >  T,  (8.11) 

for  which 

h^z)  =  exp{vT[exp(1/ao2z2)  -  1]}  (8.12) 

and  (ii)  the  sinusoidal  signal  in  (8.3),  for  which 

h2(z)  =  exp  {vT[exp(y4o2z2)I0(y4o2z2)  -  1]}.  (8.13) 

We  wish  to  compare  the  cumulative  distributions  of  v  for  these  two  signals,  and 

2  P 

we  choose  o  so  that  both  processes  v(t)  have  the  same  variance  vT,  taking  o  =  1 

2 

for  the  first  and  o  =2  for  the  second.  Thus 

z^/2 

h1 (z)  =  exp  vT(e  -  1  ) ,  (8.14) 

2 

h2(z)  =  exp  vTCeZ  /2IQ(z2/2)  -  1].  (8.15) 

Again  taking  V  <  0,  so  that  the  principal  saddlepoint  xQ  lies  on  tne 
positive  real  axis,  we  find  that  for  signal  (i)  the  path  of  steepest  descent 
passing  vertically  through  Xq  bends  around  and  approaches  zero  along  the 
hyperbola  xy  =  ir.  Indeed,  the  phase  of  the  integrand  in  (1.9)  is  now 

z2/2 

4'1  (z )  =  vT(e  -  1)  -  |  V  j  z  -  In  z,  (8.16) 

and  along  the  path  of  steepest  descent 

2  2 

Im  y^z)  =  vT  exp  V2(x  -y  )  sin  xy  -  |v|y  -0  =  0, 

z  =  x  +  iy  =  reiS.  (8.17) 

As  the  point  z  moves  to  the  right,  the  exponential  function  increases  rapidly, 
and  sin  xy  must  go  to  zero  in  order  to  keep  Im  Y^z)  equal  to  zero;  that  is, 
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xy  •»  it.  All  the  higher  branches  of  the  path  of  steepest  descent  in  the  upper 
halfplane  are  thus  bent  hairpins  asymptotic,  as  x  ®,  to  the  hyperbolas 

xy  =  (2n+1)ir,  n  =  0,  1,  2,  ...  .  (8.18) 

Each  such  branch  possesses  a  saddlepoint  of  the  integrand.  For  signal  (ii)  the 
paths  of  steepest  descent  exhibit  a  similar  configuration.  Numerical  integration 
along  such  a  complicated  set  of  curves  would  be  quite  lengthy  and  involved. 

The  real  part  of  the  phase  y 1 ( z )  is 

Re  ^(z)  =  vT  exp  14(x2-y2)  cos  xy  -  vT  -  (Vjx  -  In  r,  (8.19) 

and  the  absolute  value  exp  [Re  ^(z)]  of  the  integrand  in  (1.9)  vanishes  a3  z 
goes  off  toward  infinity  anywhere  in  the  sector 

tt/4  <  arg  z  <  v/2. 


The  originally  vertical  contour  of  integration  in  (1.9),  with  c  =  xQ  >  0,  can 
therefore  be  deformed  into  a  contour  passing  through  the  saddlepoint  xQ  and 
approaching  the  line  x  =  y  asymptotically.  For  this  reason  we  have  chosen  to 
integrate  (1.9)  along  a  hyperbolic  contour  on  which 

z  =  (xQ2  +  y2)/z  +  iy,  (8.20) 


whereupon  (1.9)  becomes 

q"(V) 


-1 

7T 


z*dy /x . 


(8.21) 


The  trapezoidal  rule  was  applied  to  integrating  this  numerically,  and  for  the 
values  of  vT  considered,  the  trapezoidal  sum  converged  rapidly  enough  for  the 
numerical  integration  to  be  feasible. 

In  Figs.  8.1  and  8.2  we  plot  the  complementary  cumulative  distributions 
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i 

q+(V)  for  these  two  signals  for  V  >  0  and  values  of  vT  ranging  from  1  to  50.  The 
I  -V 

I  abscissa  represents  the  variable  (vT)  2V.  The  ordinates  were  converted  by  the 

transformation  -erfc_,(q+),  under  which  a  normal  distribution  plots  as  a  straight 
line.  The  dashed  line  marked  '  ■»'  indicates  the  asymptote  that  the  curves 
approach  as  with  increasing  vT  the  distribution  of  the  clutter  amplitude  becomes 
more  nearly  Gaussian.  For  the  lowpass  signal  (i)  the  approach  to  normality  is 
more  rapid  than  for  the  sinusoidal  rectangular  pulse  (ii). 
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Complementary  cumulative  distribution  of  clutter  noise  consistin 
narrowband  rectangular  amplitude  modulations  of  a  carrier  having 
Gaussian-distributed  amplitudes.  Curves  are  indexed  with  the  va 

vT. 


II.  Integer-Valued  Random  Variables 
9.  Calculating  Distributions  from  Probability  Generating  Functions 
Let  w  be  a  positive  integer-valued  random  variable  with  a  probability 


distribution 


Pr  (w=k)  =Pki  k=0,  1,  2,  ...» 


pk  =  11 


with  cumulative  distribution  function  (c.d.f.) 


Q_  =  Pr  ( 


w  <  n)  V 


(9.1  ) 


and  with  complementary  cumulative  distribution  function  (c. c.d.f.) 


Qn  -  1  '  Qn  '  Pr  (w  ^  n)  =2-r  P k’ 


(9.2) 


In  typical  applications  w  is  the  number  of  electrons  counted  during  a  specific 
interval  (0,  T)  at  the  output  of  some  photoelectric  device.  The  c.d.f.  and  the 
c. c.d.f.  can  be  computed  from  the  probability  generating  function  (p.g.f.) 


h(z )  =  E(z 


(9.3) 


by  the  contour  integrals 


,  -  /  z  n  h(z)  dz 

n  =  /  1  -  z  2iti  * 

*/C 

+  /  z  n  h(z)  dz 

n  J  z  -  1  2ni  ’ 


(9.4) 


(9.5) 


'  v.\  -- 


in  which  C_  and  C+  are  circles  centered  at  the  origin  z  =  0;  C+  also  encloses  the 
point  z  =  1,  C_  does  not;  and  the  contours  enclose  no  singularities  of  the  p.g.f. 
h(z),  all  of  which  lie  outside  the  unit  circle. 

As  shown  in  [333,  the  numerical  evaluation  of  these  integrals  is  most  effi¬ 
cient  when  the  contours  pass  through  saddlepoints  of  the  integrand,  that  is, 
through  points  Xq  that  are  solutions  of 


where 


Y'  (z) 


h ’ (z )  n  1 

h(z)  z  z  -  1 


Y(z)  =  In  h(z)  -  n  In  z  -  ln[±(z  -  1  )] 


(9.7) 


is  the  "phase"  of  the  integrand,  ±(z  -  1)  standing  for  1  -  z  when  (9.9)  is  inte¬ 
grated  and  for  z  -  1  when  (9.5)  is  integrated.  One  such  saddlepoint,  xQ  =  Xg, 
lies  in  0  <  xQ  <  1,  the  other  lies  on  the  real  z-axis  between  x  =  1  and  the  left¬ 
most  singularity  of  h(z).  The  saddlepoint  can  usually  be  most  expeditiously 
located  by  solving  (9.6)  by  Newton's  method  or  the  secant  method.  Alternatively, 
the  equivalent  of  (1.17)  can  be  used,  or  y(x)  can  be  minimized  in  the  appropriate 
range  of  x  =  Re  z. 

The  circular  contours  in  (9.9)  and  (9.5)  often  suffer  the  disadvantage  that 
the  integrand  is  oscillatory  along  them,  so  that  too  many  steps  may  be  needed  in 
order  to  obtain  an  accurate  value  of  Qn+  or  Qn~  by  numerical  integration.  The 
value  of  the  integrand  decreases  most  rapidly  from  its  value  at  the  saddlepoint 
when  the  contour  is  taken  as  the  path  of  steepest  descent  through  the  saddle¬ 
point.  Along  that  path  Im  Y(z)  =  0.  Determining  the  path  of  steepest  descent, 
however,  would  much  protract  the  evaluation  of  the  contour  integral.  In  many 
cases  it  has  been  found  that  the  path  of  steepest  descent,  which  lies  symmetri¬ 
cally  about  the  real  axis,  has  nearly  a  parabolic  form,  and  it  is  then  expedi¬ 


tious  to  use  instead  its  osculatory  parabola,  that  is,  a  parabola  passing  through 


the  saddlepoint  zQ  and  having  the  same  curvature 


k  =  3  »"’(20)/»"(z0) 


(9.8: 


at  the  point  z  =  zQ.  (Primes  indicate  differentiation. )  Along  this  parabola 


z  =  zQ  +  1/2  <y  +  iy, 


(9.9) 


and  the  integral  to  be  evaluated  numerically  is 


/ 

D  h(z)z  n  . 

Re  Hz  -  1)  0  '  licy)  dy 


(9.10) 


The  initial  step  size  is  taken  as 


Ay  = 


[v  (zQ)] 


-1  /2 


(9.11  ) 


and  successively  halved  until  the  values  of  the  probability  stabilize  to  the 
number  of  significant  figures  desired.  The  trapezoidal  rule  is  used  for  reasons 
described  by  Schwartz  [6]  and  Rice  [31. 

In  other  problems,  such  as  the  one  to  be  treated  in  Sec.  10,  the  path  of 
steepest  descent  may  be  more  complicated.  It  is  advisable  first  of  all  to  trace 
a  few  such  paths  for  typical  parameter  values  by  the  method  of  Sec.  1(b).  If  one 
knows  no  simple  alternative  path  along  which  the  integral  converges  sufficiently 
rapidly,  one  can  integrate  along  a  close  approximation  of  the  path  of  steepest 
descent  by,  for  instance,  utilizing  the  five-point  quadrature  formula  (1.27)  as 
described  in  Sec.  1(d). 
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10.  Single-Stage  Photomultiplication 

The  contour-integration  method  has  been  applied  to  finding  the  distribution 
of  the  number  of  electrons  at  the  output  of  a  photomultiplier  with  a  single  stage 
of  multiplication  [34],  Light  striking  a  photoelectrically  emissive  surface 
ejects  a  number  k  of  photoelectrons  during  an  interval  (0,  T)  with  probability 
nk .  The  p.g.f.  of  this  distribution  is 


00 


k=0 


Each  such  "primary"  electron  is  accelerated  by  an  applied  voltage,  impinges  on  a 

second  metallic  surface,  and  ejects  from  it  a  random  number  of  secondary 

( s ) 

electrons.  Let  p^  be  the  probability  that  a  given  primary  electron  ejects  m 
secondary  electrons;  the  p.g.f.  of  this  distribution  is 


g(z) 


00 

E 


n  (s)  111 
p  z  . 
Km 


(10.2) 


Then  the  p.g.f.  of  the  total  number  of  secondaries  during  (0,  T)  is 


h(z)  =  pnzn  =  f  (g(z) ) , 
n=0 


(10.3) 


where  pn  is  the  probability  that  n  secondary  electrons  appear  at  the  output  of 
this  device. 

In  this  work  it  was  assumed  that  the  distribution  of  the  number  of  secondary 
electrons  per  primary  has  the  Poisson  form,  so  that 

g(z)=e^z1\  (10. U) 


where  G  is  the  gain.  The  output  distribution  can  then  be  computed  from 


(10.5) 


Pn=2^nk(kC)  e  /n ! 
k=1 

with  the  cumulative  distribution  determined  as  in  (9.1)  by  addition. 
Alternatively,  one  can  use  the  finite  relations 

pQ  =  f (e  G) , 
n 

Pn  =  Gn^  fk  S(n,  k)  e~kG,  n  >  0, 
k  =  1 


(10.6) 


where 


f 


k 


(10.7) 


and  where  S(n,  k)  are  modified  Stirling  numbers  of  the  second  kind  and  obey  the 
recurrent  relations 


SO  ,  1  )  =  1 ,  S(k,  n)  =  0,  k  >  n, 

(10.8) 

S(n  +  1,  k)  =  [S(n,  k  -  1)  +  kS(n,  k)j/(n  +  1). 

These  formulas  were  used  to  compute  exact  values  of  the  probabilities  and  Q 

for  comparison  with  those  computed  by  our  approximation  methods.  They  involve 

lengthy  computations  when  the  number  n  is  large,  and  round-off  error  introduces 

+  ~  — 

large  relative  error  into  the  c.c.d.f.  =  1  -  Qn  when  Qn  is  computed  as  in 
(9.D. 

Three  types  of  primary  distribution  j n k }  were  treated:  (i)  that  arising 
when  the  incident  light  has  a  Lorentz  spectral  density,  and  for  which  the  proba 
bilities  n  can  be  calculated  by  a  method  given  by  Bedard  [35J;  (ii)  the  nega¬ 


tive  binomial  distribution 


v  =  N  /(N  +  M)  , 

P  P 


(10.9) 


n 


k 


r(M  1  k)  , 
k'.r(M)  U 


k 

v  , 


with  Np  the  mean  number  of  primary  electrons  and  M  the  number  of  temporal  modes 
or  "degrees  of  freedom",  given  roughly  by  the  product  of  the  bandwidth  W  of  the 
incident  light  and  the  duration  T  of  the  counting  interval;  and  (iii)  the  Poisson 
distribution 


n ,  =  N  k  exp  (-N  )/k! ,  (10.10) 

k  p  p 

for  which  the  output  distribution  possesses  the  Neyman  Type-A  form. 

For  primary  distributions  of  type  (i)  and,  for  M  an  integer,  of  type  (ii), 
the  p.g.f.  f(z)  possesses  poles  z^  lying  on  the  positive  real  z-axis  to  the  right 
of  z  =  1 ,  whereupon  the  output  p.g.f.  h(z)  possesses  poles  at  the  points 

=  1  +  G_1(ln  zk  +  2rmi  ) ,  (10.11) 


for  all  positive  and  negative  integral  values  of  r,  including  r  =  0.  The  contour 
integral  (9.5)  can  then  be  evaluated  as  a  residue  series  with  terms  corresponding 
to  each  point  of  this  doubly  infinite  lattice  of  poles.  For  the  negative-binomi¬ 
al  primary  distribution,  the  array  shrinks  to  a  vertical  column  of  multiple  poles 
of  order  M  at  the  points 

C  =  1  +  G  (-In  v  +  2riri).  (10.12) 


When  the  number  n  of  output  electrons  is  of  the  order  of  or  greater  than  the 
expected  value 


E(n)  =  N  G,  (10.13) 

only  a  very  few  of  the  poles  above  and  below  the  real  axis  contribute  signifi¬ 
cantly  to  the  residue  series.  The  number  of  columns  of  the  lattice  that  need  to 
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be  included  under  case  (i)  (incident  light  with  a  Lorentz  spectral  density)  is  cf 
the  order  of  the  time-bandwidth  product  WT.  Formulas  for  computing  the  residues 
and  numerical  examples  are  given  in  [34].  Unless  WT  is  very  large,  this  residue 
method  is  most  efficient  for  computing  the  c.c.d.  Q^+  with  n  >  E(n). 

Consideration  was  given  to  using  the  contour  integrals  in  (9.4)  and  (9.5) 
for  computing  the  c.d.  and  the  c.c.d.  of  the  number  of  output  electrons.  The 
path  of  steepest  descent  is  now  too  complicated  to  be  utilized  or  even  approxima¬ 
ted,  for  it  consists  of  a  number  of  hairpin-like  branches  opening  to  the  right 
and  passing  around  each  of  the  poles  of  h(z).  We  therefore  used  instead  a 
straight  vertical  path  of  integration  passing  through  the  saddlepoint  of  the 
integrand  lying  on  the  real  z-axis;  the  one  lying  in  0  <  z  <  1  was  used  for 
,  and  the  one  lying  in  1  <  z  <  5^°^  was  used  for  Qn+>  these  having  been 
determined  by  solving  (9.6)  by  Newton’s  method.  A  bound  on  the  error  incurred  by 

truncating  the  numerical  integration  was  worked  out.  Values  of  Q  and  Q  +  over 

n  n 

a  broad  range  of  values  of  n  could  be  accurately  computed  in  this  way  without 
requiring  more  than  one  hundred  or  so  steps  of  numerical  integration.  Details 
are  given  in  Sec.  4  of  [34].  This  method  does  not  require  that  the  number  M  in 
(10.9)  be  an  integer.  The  contour-integration  method  can  be  used  for  all  three 
types  of  primary  distribution. 

On  each  branch  of  the  path  of  steepest  descent  of  the  integrand  in  (9.4)  and 
(9.5)  there  is  a  point  zk  at  which  the  magnitude  of  the  integrand  is  maximum; 
this  is  a  saddlepoint  of  the  integrand.  The  contribution  of  each  branch  to  an 
integration  of  (9.4)  or  (9.5)  along  its  path  of  steepest  descent  can  be  crudely 
approximated  by  assuming  that  in  the  neighborhood  of  each  saddlepoint  zk  the 
integrand  has  a  Gaussian  form,  and  this  leads  to  the  formula 

00 

Qn±  *  ^  ^  [2n'P"(zk)]  1/2  exp  [y(zk)],  (10.14) 

k=-» 
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where  'y(zj<)  is  the  phase  of  the  integrand  as  in  (9.7)  and  Y"(zk)  is  its  second 
derivative  evaluated  at  the  k-th  saddlepoint  [3 j -  For  the  most  part  two  or  three 
saddlepoints  on  each  side  of  the  principal  one  xQ  on  the  real  z-axis  are  all  that 
it  is  worth  taking  into  account;  often  the  one  at  xQ  alone  suffices.  Improved 
accuracy  was  obtained  by  evaluating  the  contribution  of  the  principal  branch  of 
the  path  of  steepest  descent,  which  crosses  the  real  axis  at  zQ,  by  means  of  the 
uniform  asymptotic  expansion  described  in  [36]  and  [37].  The  details  of  tms 
method  are  to  be  found  in  Sec.  5  of  [3^3- 


11.  Multistage  Photomultiplier 


A  photomultiplier  usually  has  not  one,  but  several  stages  of  electron  multi 

plication.  We  number  the  stages,  each  associated  with  an  electrode  called  a 

"dynode,"  from  the  final  one.  A  single  electron  striking  the  k-th  dynode  ejects 

( k) 

n  secondary  electrons  from  it  with  probability  p^  ;  the  collection  of  these 
probabilities  constitutes  the  k-th  "single-stage"  distribution  and  has  a  p.g.f. 

00 

gk(z)  =  pnk)  zn*  (1K1) 

n=0 

Then  if  there  are  N  stages  in  all,  the  probabilities  pn  that  n  electrons  emerge 
at  the  output  when  a  single  primary  electron  strikes  the  N-th  dynode  have  the 
p.g.f. 


where 


H»U) 


00 


m=0 


(1  1 .2) 


H1 (z)  =  g1 (z) , 

Hk(z)  =  gkK-1(z)^  k  =  2-  3 . N. 


(11 .3) 


If  the  distribution  of  primary  electrons  again  has  the  p.g.f.  f(z)  defined  in 
(10.1),  the  total  number  of  output  electrons  has  the  p.g.f.  f(H^(z)). 

In  this  study  we  have  concentrated  on  calculating  the  cumulative  output 
distributions,  defined  as  in  (9.1)  and  (9.2),  for  a  single  primary  electron,  so 
that  in  (9.3)-(9.5)  the  p.g.f.  of  the  output  is 

h (z )  =  H  (z  ) . 

N 

We  assumed  that  all  the  single-stage  distributions  have  the  same  negati ve-binomi 


al  or  Polya  form,  with  common  p.g.f. 

gk(z)  =  g(z)  =  [1  -  bG(z  -  l)]  1/b,  b  >  0,  (11.4) 

which  corresponds  to  the  distribution  in  (10.9)  with 

M  =  1/b,  v  =  bG/( 1  +  bG). 

Again  G  is  the  gain  per  stage;  the  overall  gain  of  the  photomultiplier  is 

Gq  =  GN,  (1 1  .5) 

where  N  is  the  number  of  stages.  The  Poisson  single-stage  distribution  represen¬ 
ted  by  (10.4)  arises  in  the  limit  b  ■*  0. 

The  singularities  of  the  overall  p.g.f.  H^(z)  now  lie  on  or  close  to  the 
portion  of  the  real  z-axis  lying  to  the  right  of  the  point  z  =  1,  and  the  path  of 
steepest  descent  has  roughly  the  form  of  a  parabola  lying  symmetrically  about  the 
real  z-axis  and  opening  to  the  right.  This  path  can  be  replaced  by  its  oscula- 
tory  parabola,  whose  curvature  <  at  the  saddlepoint  is  determined  by  (9.8).  The 
first  three  derivatives  of  the  phase  Y(z)  of  the  integrand,  defined  as  in  (9.7), 
are  computed  by  an  N-fold  set  of  three  recurrent  relations.  The  cumulative  prob¬ 
abilities  Qn+  are  then  computed  by  numerical  integration  of  (9.10).  The  results 
were  compared  with  values  computed  by  the  recurrent  relations  given  by  Prescott 
[38]  for  N  «  3  and  values  of  n  up  to  100,  and  the  agreement  was  excellent  [39]. 

The  recurrent  relations  require  of  the  order  of  n  operations  to  compute  pn  or 
+ 

Qn  ,  and  this  number  is  independent  of  the  number  of  significant  figures 
sought.  Thus  for  n  *  3000  of  the  order  of  nine  million  operations  would  be 
required,  in  contrast  to  the  ten  or  twenty  steps  of  numerical  integration  of 
(9.10)  needed  for  adequate  precision. 

+ 

When  the  cumulative  distributions  Qn  so  calculated  are  plotted  on  a  semi- 
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logarithmic  grid  for  various  numbers  N  of  stages  of  multiplication,  the  overall 
gain  Gq  =  GN  being  kept  fixed,  one  finds  the  points  lying  closer  and  closer  to 
straight  lines  as  N  increases.  The  results  are  fitted  quite  closely  by  a 
geometric  or  "exponential"  distribution, 

pQ  =  1  -  (1  -  v)G0, 

_  * .  x2  n-1  .  « 

pn  -  G0O  -  *)  v  ,  n  >  0,  („.6) 

Q„‘  -  G0(1  -  v)*"-',  n  >  0, 
in  which  the  parameter  v  is  given  by 

k2(G  -  1) 

v  =  a/(a  +  1),  a  =  gGfG"-7  i )  ~*  (11.7) 

where  <2  is  the  second  factorial  moment  of  the  single-stage  distribution  and  G  is 
the  gain  per  stage, 

G  =  g’(1),  «c2  =  g"(1),  (11.8) 


primes  denoting  differentiation.  This  value  of  v  is  chosen  so  that  the  distribu¬ 
tion  in  (11.6)  has  the  same  variance  as  the  exact  distribution, 


Var  n 


°VGo  -  1 ) 

G(G  -  1) 


(11.9) 


A  number  of  observers  have  reported  approximately  exponential  distributions 
for  the  number  of  output  electrons  in  photomultipliers,  but  an  adequate  explana¬ 
tion  of  this  phenomenon  seems  to  have  been  lacking  [38].  In  [39]  the  distribu¬ 
tion  in  (11.6)  was  shown  to  approach  the  true  output  distribution  under  the  con- 

N 

ditions  stated:  the  overall  gain  Gq  =  G  is  fixed,  the  single-stage  distribu¬ 
tions  are  identical,  and  the  number  N  of  stages  increases,  the  gain  G  per  stage 
decreasing  simultaneously  toward  1 .  A  related  limit  theorem  in  the  theory  of 
branching  processes  has  been  derived  by  Fahady,  Quine,  and  Vere-Jones  [HO]. 
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12.  Performance  of  Receivers  with  Avalanche  Photodiodes 
(a)  Distribution  of  the  Output  of  an  Avalanche  Photodiode 

In  an  avalanche  photodiode  primary  electrons  generated  at  one  face  by 
incident  light  are  accelerated  under  an  applied  voltage  and  create  hole-electron 
pairs  in  the  body  of  the  diode  by  collision.  The  holes  and  the  secondary 
electrons  are  themselves  accelerated  and  create  more  such  pairs.  The  probability 
P r  (n|  N)  that  a  fixed  number  N  of  primary  electrons  generate  a  total  of  n 
electrons  at  the  output  was  calculated  by  McIntyre  under  the  assumption  that  the 
collision-ionization  probability  for  holes  is  a  fixed  fraction  g  of  that  for 
electrons,  independently  of  their  energy,  0  <  g  <  1  [41].  Personick  [42]  showed 
under  the  same  assumption  that  the  probability  generating  function  (p.g.f) 

00 

M(z)  =  ^  Pr  (n |  1)  zn  (12.1) 

n=0 

of  the  distribution  of  the  number  of  output  electrons  arising  from  a  single 
primary  electron  is  the  solution  of  the  equation 

z  =  M[  1  +  a(M-1)]"b,  b  =  (1-g)_1  >  1,  (12.2) 

where  a  is  related  to  the  gain  G  of  the  diode  by 

_  5 

G  =  ( 1 -ab )  ' ,  0  <  a  <  1  .  (12.3) 

Ordinarily  the  number  of  primary  electrons  has  a  Poisson  distribution  with 
mean  X  proportional  to  the  intensity  of  the  incident  light,  and  the  number  n  of 
electrons  at  the  output  has  a  compound  Poisson  distribution, 

V*  xNe~* 

pn  *  Zw  Vr~  *  (nl  N)*  (12-4) 

N=0 
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The  cumulative  distribution,  that  is,  the  probability  that  fewer  than  k 
electrons  are  counted,  must  then  be  computed  by  summing  the  probabilities  pn, 

k-1 

n=0 

and  the  complementary  cumulative  distribution  is  qk  +  ■  1  -  q  .  Computing  these 
probabilities  in  this  manner  is  tedious  because  of  the  complexity  of  McIntyre's 
formula  [41]  for  P r  (nj  N),  and  when,  as  usually  in  a  multiplicative  device  such 
as  this,  the  numbers  k  are  large,  many  summations  are  required  and  round-off 
error  is  serious. 

The  method  of  saddlepoint  integration  has  been  found  most  efficient  for 
determining  the  cumulative  distribution  qk  of  the  number  of  output  electrons 
[43].  One  computes  qk  or  qk+  by  numerical  integration  of  (9.4)  or  (9.5),  in 
which  now  the  probability  generating  function  is 

h(z)  =  exp[A(M-1 )],  (12.6) 


M  and  z  being  related  by  (12.2).  Rather  than  solving  (12.2)  for  M  at  each  point 
z  on  the  contour  of  integration,  we  change  variables  in  (9.4)  and  (9.5)  and  use  M 
as  the  variable  of  integration,  obtaining  the  formula 


A(M-1  ) 
e 


M  n[1  ♦  a(M-1)]bn-1 


r  1  -  a  +  a ( 1  - b ) M  1  dM 

*  [ - d  5-7,  U2.7) 

M  -  [ 1  +  a(M-1 )]° 

where  C_  and  C+  are  curves  in  the  M-plane  passing  through  the  image  MQ  or  MQ+  of 
the  saddlepoint  x^  or  xQ+  of  the  integrand  in  (9.4)  or  (9.5).  The  z-plane  is 


cut  along  the  real  z-axis  from  a  point 


=  (1-a)/[a(b-1 )]  >  1, 


to  z  »  »,  in  order  to  render  the  solution  of  (12.2)  for  M  single-valued.  A 
parabolic  contour  in  the  cut  M-plane  was  found  to  yield  the  probabilities  qk+  and 
qk~  to  seven  significant  figures  with  fewer  than  twenty  steps  of  numerical 
integration,  even  for  values  of  k  as  high  as  50000.  Indeed,  the  larger  the  gain 
G  and  the  mean  number  \  of  primary  electrons,  the  more  accurate  this  method  of 
saddlepoint  integration  becomes. 


fJ 


(b)  Distribution  of  the  Sum  of  Diode  Output  and  Gaussian  Noise 

In  an  optical  receiver  utilizing  an  avalanche  photodiode,  Gaussian  noise 
corrupts  the  diode  output  in  a  subsequent  amplifier.  The  performance  of  such  a 
receiver  was  calculated  by  Personick,  Balaban,  Bobsin,  and  Kumar  [44]  by 
convolving  the  distribution  of  the  number  of  output  electrons,  as  computed  by 


(12.4),  with  a  Gaussian  distribution  representing  the  noise.  This  computation 
was  so  lengthy  that  they  resorted  instead  to  Monte  Carlo  simulation,  crude 
Gaussian  approximations,  and  Chernoff  bounds.  The  distribution  of  the  output  can 
be  computed  easily  by  the  method  of  saddlepoint  integration. 

The  moment-generating  function  (m.g.f)  of  the  sum  v  of  the  number  of  output 

2 

electrons  and  a  Gaussian  noise  variate  with  mean  zero  and  variance  o  is 


I 


I 


-V.  We  used  the  method  of  Sec.  1(d). 

One  computes  an  approximate  path  of  steepest  descent  passing  through  the 
saddlepoint  of  the  integrand  of  (1.9)  or  (1.10)  by  the  method  described  in  Sec. 
1(b).  The  integration  is  carried  out  along  this  segmented  path, 


00  /*M. 


q~(v;  =  /  sW)  dM/7rl» 

k=0  "M, 


(12.10) 


where  we  have  replaced  the  integration  variable  z  by  M,  and 


g(M)  =  ( ±s )  1  expU(M-l)  +  ,4o2s2](^)  1 


(12.11) 


with 


dM  _  M[ 1  +  a(M-1 )] 
ds  1  +  a(M-1 )  -  abM 


and  with  s  given  by  (12.9)  as  a  function  of  M.  The  saddlepoint  sQ  =  fg(M0)  is 
obtained  by  solving  the  equation 

f ' (M)  =A^+o2s-V--=0,  s=f  (M),  (12.13) 

ds  s  s 


for  M  =  Mq  by  Newton's  method.  The  contribution  of  each  term  in  (12.10)  is 
evaluated  by  the  five-point  formula  (1.27)  of  Birkhoff  and  Young  [7].  The 
initial  length  j  AM f  =  6  of  the  segments  can  be  taken  as 


6  -  %  #  ) 

2  ds 


M=M0'-'1'"(S0)]  2* 


(1  2.U) 


where  Sq  =  fg(MQ)  is  the  saddlepoint  on  the  Re  s-axis;  the  second  derivative 
¥"(s0)  will  have  been  calculated  in  determining  the  saddlepoint  by  Newton's 
method. 

One  can  also  apply  the  saddlepoint  approximation  described  in  [8]  and  [46] 
to  this  problem,  evaluating  the  tail  probabilities  q~(V)  and  q+(V)  as 
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q"(V)  =  [2ir’(',! (Sq)  ]  2  exp  ^CSq)  , 


where  the  phase  <P(s)  is,  as  in  (1-15), 


»(s)  =  X  (M-1  )  +  V2o2s2  -  Vs  -  In  (±s), 


(12.15) 


(12.16) 


and  the  saddlepoint  Sq  is  the  root  of  (12.13)  lying  to  the  right  of  the  origin 
for  q  +  (V)  and  to  the  left  for  q“(V).  By  using  M  as  the  independent  variable, 
(12.13)  can  be  solved  by  Newton's  method,  in  the  course  of  which  i"'(sq)  is 
determined . 

In  the  binary  receiver  analyzed  by  Personick  et  al .  [44]  tne  decision  level 
V  above  which  it  decides  that  a  signal  is  present  is  set  to  minimize  the  average 
probability  of  error.  As  0's  and  1 's  were  taken  as  equally  frequent,  this  is  the 
point  at  which 


Pq(V)  =  p,(V), 


(12.17) 


where  Pq(v)  and  P-|(v)  are  the  probability  density  functions  of  the  output  v  under 
the  hypotheses  H0  (no  signal)  and  H.]  (signal  present),  respectively.  The  mean 
numbers  A  ^  of  input  electrons  under  these  two  hypotheses  are 


X1  =  Ns  +  V 
X0  =  +  V 


(12.18) 


3 


where  Nd  measures  the  strength  of  the  dark  current,  Ng  the  strength  of  the  signal 
under  H.|  ,  and  n  <<  1  the  leakage  in  the  optical  transmitter.  The  aim  of  the 
computations  described  by  Personick  et  al .  [44]  was  to  determine  the  minimum 
Input  signal  strength  Ng  and  the  associated  decision  level  V  so  that  the  average 
probability  of  error 


p  =  l4tqn  (V)  +  q  "(V)] 


(12.19) 


takes  on  a  pre-assigned  value.  Here  the  subscripts  on  A,  p,  and  q-  refer  to  the 


hypothesis. 

The  density  functions  Pq(V)  and  p^V)  needed  for  setting  the  decision  level 
V  through  (12.17)  can  also  be  calculated  approximately  by  the  saddlepoint  method 
[8,  46], 

p  (V)  =  h  (a1)  exp  (-s.V)  [2u<J>."(si)]”1/2,  i  =  0,  1,  (12.20) 

where  h^s)  is  the  m.g.f.  of  the  output  v  under  hypothesis  H.  and 

4>  (s)  =  In  h  ( s)  -  sV  =  A.(M-I)  +  ‘4o2s2  -  sV  (12.21) 

is  the  associated  modified  phase.  The  saddlepoints  s .  are  the  roots  of  the 
equation 

♦  .'(s)  =  A.  ^  +  o2s  -  V  =  0,  i  =  0,  1,  (12.22) 

l  l  ds 

which  can  be  solved  by  a  few  iterations  of  Newton's  method  for  the  saddlepoints 
sQ  and  s1 ,  again  treating  M  as  an  independent  variable  and  determining  s  from 
(12.9).  The  equation  (12.17)  for  the  decision  level  V  is  in  this  approximation 
equivalent  to 

A  i  ( H 1  —  1 )  ♦  V2o2s  1 2  -  s  V  -  l41n  <t  1 "  ( s  1 ) 

A0(Mo"1)  +  Vao  sQ  -  sQV  -  41n  $0"(s0), 

s.  =  f  (M. ) ,  i  =  0,  1 .  (12.23) 

1  si 

This  equation  is  solved  for  V  by  the  secant  method. 

One  then  uses  this  value  of  the  decision  level  V  in  (12.19),  computing 
qg+(V)  and  q1*(V)  by  the  saddlepoint  approximation  (12.15)  in  order  to  obtain  the 
error  probability  as  a  function  Pe(N„)  of  the  input  signal  strength  N  .  By  using 


the  secant  method  one  can  then  determine  the  signal  strength  Ng  to  yield  the  pre¬ 
assigned  value  of  the  error  probability  Pe. 

The  optimization  procedure  described  in  [^5]  starts  with  Gaussian  approxi¬ 
mations  to  the  densities  p^(V)  and  the  cumulative  probabilities  Qq+(V)  and 
q -I - ( V ) .  The  means  and  variances  of  v  needed  for  these  are  easily  calculated,  and 
the  equation  (12.17)  for  the  decision  level  V  reduces  in  this  approximation  to  a 
quadratic  equation.  One  uses  the  secant  method  to  solve  (12.19)  in  this  Gaussian 
approximation  to  obtain  starting  values  of  Ng  and  V  for  the  computation  with  the 
saddlepoint  approximations  just  outlined.  Two  or  three  iterations  of  each 
suffice  to  produce  accurate  estimates  of  the  minimum  signal  strength  Ng  and  the 
decision  level  V  needed  to  attain  the  pre-assigned  false-alarm  probability  P  . 

In  a  further  stage  of  this  optimization  procedure  the  probabilities  qQ+(V) 
and  q-|_(V)  were  computed  by  the  numerical  saddlepoint  integration  technique 
described  earlier  in  this  section.  The  decision  level  was  still  calculated  by 
the  saddlepoint  approximation  (12.20)  -  (12.23).  Although  one  could  also 
calculate  the  density  functions  Pq(v)  and  P^(v)  more  accurately  by  numerical 
integration  of  (1.7),  this  refinement  produced  so  slight  a  change  in  the  final 
results  as  not  to  be  worth  the  trouble.  In  Table  12.1  we  list  the  values  of  the 
signal  strength  Ng  and  the  decision  level  V  determined  by  the  three  successive 
parts  of  this  optimization  technique.  It  is  seen  that  the  third  part,  requiring 
numerical  quadrature,  alters  the  values  of  Ng  and  V  very  little  from  those 
yielded  by  the  saddlepoint  approximations  in  the  second  part.  For  most  practical 
purposes  carrying  the  procedure  through  only  the  first  and  second  parts  suffices 
and  yields  results  far  more  accurate  than  those  of  the  Gaussian  approximation 
alone,  but  without  much  additional  computation. 
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13.  Photoelectron  Counting 

(a)  Incident  Light  of  Arbitrary  Spectral  Density 

When  light  falls  on  a  photoelectr ically  emissive  surface,  the  probability 
generating  function  of  the  number  of  electrons  ejected  during  an  interval  (0,  T) 
is  given  by 


h(z)  =  E[eI(z_1 }], 


(13-1  ) 


where  I  is  a  random  variable  proportional  to  the  total  energy  falling  on  the 
surface  during  that  interval.  If  we  assume  the  light  to  be  plane-polarized, 
quasimonochromatic ,  and  spatially  coherent  over  the  entire  emissive  surface,  its 
field  can  be  written 


v ( t )  =  Re  [V(t)e  w  ],  0  <  t  <  T, 


in  which  the  complex  envelope  V(t)  is  a  circular  complex  Gaussian  random  process, 
and  n  =  2irv  is  the  central  angular  frequency  of  the  light.  Then 


•'a 


t)  I  dt 


(13.2) 


with  n  a  constant  directly  proportional  to  the  quantum  efficiency  of  the  emissive 
surface  and  inversely  proportional  to  the  photon  energy  hv. 

The  expected  value  E[V(t)3  =  S(t)  represents  the  complex  envelope  of  a 
coherent  light  signal,  such  as  the  output  of  a  pulsed  laser.  The  difference 
( t )  =  V(t)  -  S(t)  represents  incoherent  light  such  as  is  produced  by  an 
ordinary  incandescent  source.  It  is  assumed  to  be  a  stationary  Gaussian  random 
process  with  a  temporal  coherence  function  defined  by 


G(VV  -  2  E[VVVb(t2)]- 


(13. 3; 
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The  spectral  density  of  tne  light  is 


/ 


fU)=  /  G(t)elut  dt. 


(13. 4) 


The  mean  number  nQ  of  photoelectrons  ejected  by  this  incoherent  component  is 


n-  =  n(i (0  )  T : 
u  •  - 


n  i  .*,) 


the  mean  number  n  ejected  by  the  coherent  component  is 


I  S( t ) I  dt. 


(13.6) 


The  p.g.f.  h(z)  can  be  expressed  as  [25] 


rTrT 

h(z)  =  [DU)]  '  exp{n  (z-1)[l  -  f  J  c*(t)iJi(t,uU;T)o(u)  dtdu]{, 

<  Jq 


5  =  n0(1-z), 


(13-7) 


where 


D(o  =]~]  n+Ak5) 

k  =  1 


(13.8) 


is  the  Fredholm  determinant  associated  with  the  integral  equation 


Vk(t) 


=  T 


■f 


Y(t-s)  fk(s)  ds,  Y(t) 


=  G(t)/G(0),  (13.9; 


whose  eigenvalues  are  A^ .  Furthermore,  a ( t )  is  a  normalized  signal  envelope, 


o(t)  =  S(t) 


f  |S(, 


)i2  du. 


(13-10) 
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and  \Ji(t,u;C;T)  is  the  resolvent  kernel  defined  by  the  integral  equation 


<J/(t 


,U;C;t)  +  |  j  4>(  t ,  V  ;  5. ;  t 


)Y(v-u)  dv 


-1 


=  £T  Y(t-u) ,  0  <  (t,u)  <  t  <  T. 


(13.11) 


In  the  past  the  distribution  of  the  number  of  photoelectrons  has  been 
evaluated  by  decomposing  the  field  V(t)  into  temporal  modes  by  a  Karhunen-Loeve 
expansion,  whereupon  the  p.g.f.  can  be  written  as  a  product  of  factors  involving 
the  eigenvalues  of  (13-9)  and  the  coefficients 


S 


k 


T 

f*(t)  S(t)  dt 


of  the  coherent  component.  The  distribution  corresponding  to  each  factor  or  mode 
has  the  Laguerre  form,  and  the  distribution  of  the  total  number  V  of  photo¬ 
electron  counts  is  the  convolution  of  as  many  of  these  as  necessary  [47].  Other 
modal  decompositions  can  also  be  used  [48].  (Further  references  are  to  be  found 
in  ref.  [25].)  Such  methods  require  much  computation  time  and  storage, 
particularly  when  the  product  WT  of  the  counting  time  T  and  the  bandwidth  W  of 
r(w)  is  large.  The  method  of  saddlepoint  integration  provides  much  greater 
efficiency  in  calculating  the  counting  distribution. 

When  the  light  has  a  Lorentz  spectral  density, 


2pG(0) 


rU)  =  <  y(t)  =  e 

ta  +p 


-ml 


(13.12) 


the  Fredholm  determinant  has  the  well-known  form  [24] 

D(0  =  e'm[(g  +  1)2emg  -  (g-1)Vmg]/(4g)  , 


1/2 


g  =  (1+25/m)  ,  5  =  n  (1-z),  m  =  pT. 


(13-13) 


If  the  signal  is  a  sinusoid  of  constant  amplitude  and  with  a  frequency  displaced 
from  ft  by  A,  as  we  have  shown,  the  p.g.f.  h(z)  is  [25] 


h(z) 


[DU)]  1  exp{n  (z-1)[-  _ 

g  +6 


1+6 


4g{(g+1(g-6  )e 


(g-1 )(g+62)e  mg 


2k [  (1-6  )cos  m6  - 


26  sin  m6 ] }  j  j 


2,  2  2 . 2 r  /  ,  s2  mg  ,  '  .2  -mg-, 

a  (g  +6  )  [(g+1)  e  -  (g-1)  e  ®] 


03.1*0 


with  6  =  A/u  and  the  rest  of  the  notation  as  in  (13.13).  After  putting  this  into 
the  integrals  (9.4)  and  (9.5)  we  have  computed  the  complementary  cumulative 
distribution  qk+  by  evaluating  the  integrals  numerically  along  a  parabolic 
contour  passing  through  the  appropriate  saddlepoint.  Even  for  values  of  k  of  the 
order  of  thousands,  accuracy  of  seven  significant  figures  was  attained  by  fewer 
than  twenty  steps  of  numerical  integration. 

For  more  complicated  spectra  r(w)  the  p.g.f.  h(z)  cannot  be  written  down  in 
closed  form.  It  can  be  obtained  numerically  at  points  on  the  contour  of 
integration,  however,  from  the  alternative  expression  [25] 


h(z)  =  [DU)] 


-1 


r  - 

x  exp{ns(z-1)  /  [o*( t U*)~o*(t  )][ o(t  U)-o(t  )]dx }  (13.15) 


where 


't;0  =  /  (  t  » u ;  i 

Vi 


i(t;0  =  /  iJ)(T,u,U)o(u)du 
■{) 

and  \ji(T;uU)  =  i|i(t,uU;t)  is  the  solution  of  the  integral  equation 


(13.16) 
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03.17) 


* 

\ 

4 


<Kt»u;£)  +  £  T 


V; 


£)Y(v-u)dv  =  £T  Y(t-u) , 


°  <  u  <  t  <  T. 

In  terms  of  this  function  the  Fredholm  determinant  is  given  by  [24] 


.  r\ 

exp[  /  i|i(t  ,  t ;  £)dt  J , 
J0 


(13-18) 


The  function  \Ji(t,u;£)  is  the  kernel  of  the  minimum-mean-square-error  causal 

y 

estimator  of  the  process  (£/T)  2Vfe(t)  in  the  presence  of  circular-complex  white 
Gaussian  noise  with  unit  spectral  density  provided  £  is  real  and  positive.  As 
such  it  can  be  computed  by  solving  Kalman-Bucy  differential  equations  as  in 
(5.18)  -  (5.23)  whenever  one  can  establish  a  state-space  model  for  generating 
Vb(t)  from  white  noise.  Indeed,  the  procedures  for  determining  the  distribution 
of  the  average  power  of  a  Gaussian  random  process,  described  in  Sec.  6,  can  be 
applied  to  evaluating  the  terms  in  h(z)  in  (13.15).  In  particular,  the  Fredholm 
determinant  D(£)  can  be  expressed  by  the  quotient  of  determinants  in  (6.9)  when 
the  spectral  density  r(w)  in  (13. 4)  is  a  rational  function  of  the  frequency,  like 
4>(<o)  in  (6.5).  Equation  (6.8)  for  the  Bk's  appearing  in  (6.9)  is  now 


1  +  £T  f(-ip)  =0,  £  »  n  ( 1-z), 


(13-19) 


with  r(u>)  =  rU)/G(0)  the  Fourier  transform  of  Y(t)  defined  in  (13.9).  The  form 
in  (6.9)  enables  the  integrands  of  (9.4)  and  (9.5)  for  the  cumulative  distri¬ 
bution  of  the  number  of  photoelectron  counts  and  its  complement  to  be  computed 
easily  at  the  points  on  the  contour  of  integration  in  the  complex  plane.  In  part 
(b)  of  this  section,  this  method  is  applied  to  light  with  a  spectral  density 
composed  of  two  overlapping  Lorentz  spectra,  as  in  (6.5)  with  n  =  2. 

When  the  product  WT  of  the  counting  time  T  and  the  bandwidth  W  of  the 


■  *  M.  Vk  '  M  A  »  - 


incoherent  light  is  large,  one  can  utilize  the  so-called  "Toeplitz  approximation" 
to  write  the  p.g.f.  as  [25] 


h(z)  a  expjn  ( z - 1 ) 
s 


/CO  » 00 

— —  p-  -T  /  ln[  1  +5T_1  r(oj)  ]dw/2ir} , 

uct_1tU)  211  J 


r(uj)  =  r(w)/G(o),  z  (o>) 


■/' 


o(t)e“i“tdt,  5  =  nQ( 1 -z) .  (13.20) 


This  can  be  substituted  into  (90)  or  (9.5)  for  numerical  integration  to 
determine  the  cumulative  distribution  of  the  number  of  photoelectrons 
approximately.  For  incident  light  with  a  Lorentz  spectral  density  the  result  has 
been  found  quite  close  to  the  true  distribution  yielded  by  (13-1*0  for  a  time- 
bandwidth  product  m  =  yT  as  low  as  [25]. 

(b)  Bimodal  Spectral  Density  (V.  Staggs) 

(i)  Introduction 

In  the  detection  of  incoherent  narrowband  light  by  counting  photoelectrons, 
the  counting  statistics  are  specified  completely  by  the  Frednolm  determinant,  if 
one  has  a  way  of  inverting  the  probability  generating  function  for  the  number  of 
photoelectrons.  From  [25]  and  03.7),  in  the  absence  of  a  coherent  light  signal, 
the  probability  generating  function  of  the  number  of  photoelectrons  is 

h(z)  =  [D(x) ]  \  x  =  nPT(l-z), 


where  n  is  the  quantum  efficiency  divided  by  hv,  P  is  the  power  in  the  spectrum, 
and  T  is  the  counting  interval.  The  probabilities 


+ 


and 


can  be  obtained  by  a  numerical  contour  integration  involving  h(z).  For  a 
spectrum  consisting  of  a  linear  combination  of  Lorentz  spectra,  the  function  h(z) 
can  be  computed  by  using  the  matrix  methods  described  in  Section  6.  Here  we 


shall  describe  a  numerical  procedure  written  in  FORTRAN  to  find  Qk  +  and  Qk_  by 
the  matrix  method  for  Lorentz  spectra. 

This  numerical  procedure,  named  FRED,  was  compared  with  known  methods  in  the 
case  of  a  single  Lorentz  spectrum  in  order  to  test  its  accuracy.  Furthermore, 
spectra  consisting  of  a  sum  of  two  Lorentz  spectra  were  used  to  explore  cases 
where  previously  known  techniques  do  not  apply  or  are  inefficient.  In 
particular,  the  FRED  algorithm  can  compute  accurate  results  when  the  spectrum  is 
bimodal,  the  counting  time  is  short,  and  one  examines  large  values  of  k  where  Qk+ 


is  small.  This  is  a  regime  where  earlier  approximations  are  inaccurate. 

(ii)  Description  of  the  algorithm 

The  fundamental  purpose  of  the  FORTRAN  program  FRED  is  to  mechanize  the 
matrix  procedure  given  in  Section  6.  The  Fredholm-determinant  algorithm  is  then 
used  to  perform  a  numerical  contour  integration  for  the  computation  of  Qk+  or  Qk_. 
An  important  secondary. aspect  of  generating  the  FRED  program  was  to  be  alert  to 
the  numerical  conditioning  of  the  various  steps  involved,  and  to  take  appropriate 
action  to  keep  the  results  robust  against  roundoff  error.  The  following  is  an 
overview  of  the  numerical  approach. 

A.  Fredholm  Determinant 

The  spectrum  of  the  narrowband  process  n(t)  is  expressed  as 


♦  U)  =  2  }  - — - , 

(oi-i  v*)  (w+i  u  ) 
m=  i  m  m 

as  in  (6.5),  except  that  we  have  taken  g_  *  0.  Here  y  =  y'  +  iu",  or 

m  m  m  m 

equivalently , 


where  is  the  frequency  offset  of  a  Lorentz  spectrum  with  width  parameter  vr 
and  power  fm.  The  solution  of  the  integral  equation  for  the  Fredholm  determinant 
involves  finding  the  2n  roots  of  the  equation 

1  +  |  $(-ip)  =  0.  p  =  B1 ,  ...  »&2n. 

This  is  accomplished  by  calling  subroutine  FILPAR  to  fill  the  array  P  of 
coefficients  of  the  polynomial  numerator  of  the  above  equation,  and  then  using 
IMSL  subroutine  ZCPOLY  to  find  the  roots  p.  Since  the  roots  must  be  ordered  so 

that  B-| . Bn  lie  on  the  branches  of  the  root  locus  originating  at  the  poles 

pm  of  <t> ( — ip )  having  positive  real  parts,  a  sorting  procedure  is  applied  after  the 
roots  are  found. 

Using  the  quantities  8-| ,  62>  •••*  82n,  W1  ’  y2 . y2n*  and  T>  the  matrices 

§1  ,  §2’  -i  ’  ~2’  -1’  '2’  and  1  are  Sene rated,  and  the  Fredholm  determinant  is 
evaluated  using  (6.9), 

det(D2_1M*_1  )  det  §1  D1  _1 ) 

D(z)  =  - 1 1 - — i - • 

det  (C^^ 

Within  subroutine  GETFRD,  which  computes  D(z) ,  the  two  determinants  in  the 
numerator,  TERM1  and  TERM2,  and  the  denominator,  TERM3>  are  computed 
separately.  The  quotient  TERM2/TERM3  is  computed  first  since  its  terms  are  of 
the  same  order  of  magnitude,  in  order  to  avoid  overflow  or  underflow. 

B.  Parameters  of  Integration 

In  order  to  perform  the  numerical  contour  integration,  one  must  know  the 
parameters  of  the  osculatory  parabola  which  approximates  the  path  of  steepest 
descent  as  described  in  Sec.  1(c)  and  in  [5].  As  in  (9.4)  and  (9.5),  the  tail 
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‘V  /Y 


probabilities  are 


-  I  z  k  h(z)  dz 
■'k  ”  /  1  -  z  2iri 

*c 

+  _  /  z  k  h(z)  dz 
k  J  +  z  -  1  2iri 


(13.21  ) 


(13.22) 


where  C  is  a  parabola  intersecting  the  real  axis  between  the  origin  and  the 
point  (1.0,  0.0),  and  C+  similarly  intersects  the  real  axis  between  (1.0,  0.0) 
and  the  nearest  singularity  of  h(z)  greater  than  (1.0,  0.0).  One  uses  (13.21) 
when  k  <  nQ  =  nPT  and  (13.22)  when  k  >  n^. 

The  integration  procedure  starts  by  finding  the  first  eigenvalue  lying  to 
the  right  of  the  point  z  =  1 .  Subroutine  LAMBDA  performs  this  search.  It  starts 
with  a  simple  Newton  iteration  (the  secant  method)  and  then  switches  to  a  second- 
order  extrapolation  method  using  inverse  interpolation  with  relaxation.  There  is 
a  reason  for  this  elaborateness.  Determinants  are  found  in  subroutine  GETFRD  by 
calling  the  IMSL  suroutinte  LEQT1C,  which  is  also  used  to  solve  systems  of  linear 
equations  involving  complex  matrices.  LEQT1C  always  performs  an  LU  factorization 
of  the  coefficient  matrix,  and  whenever  the  matrix  is  nearly  singular,  the 
factorization  algorithm  may  attempt  to  divide  by  a  very  small  pivot,  sometimes 
resulting  in  a  fatal  error.  Thus,  it  makes  sense  to  find  the  first  eigenvalue 
without  actually  evaluating  the  Fredholm  determinant  with  that  eigenvalue,  which 
implies  that  one  needs  a  powerful  predictor  such  as  the  second-order  technique. 
Furthermore,  in  order  to  keep  the  predictor  from  jumping  immediately  to  the  first 
eigenvalue  A,  a  relaxation  is  performed,  whereby  the  next  value  of  A  generated  is 
only  part  way  between  the  predicted  value  and  the  previous  value.  In  this  way,  a 
sequence  of  approximations  L  is  generated  that  stays  to  the  left  of  A,  and 
reasonably  far  from  it,  but  which  converges  to  A  in  a  way  that  permits  one  to 
generate  a  predicted  A  that  converges  well.  This  avoids  numerical  problems  in 


to C'  /c*:  ;s  v  / 


/  .j.  ••  ■  >  ■  -  v  -.'j.  • .  .  •  .V.\.  -j  .  • 


LEQT1C,  it  guarantees  that  X  is  approached  from  the  left  and  that  one  does  not 


jump  over  it  because  of  roundoff  error  and  erroneously  deliver  the  next 


eigenvalue,  and  it  converges  adequately  fast. 

Once  X  is  found,  one  can  compute  the  quantities  xQ,  <,  and  T'(Xq),  where  as 
in  (9.7) 


'i'(z)  =  In  [±z  kh(z)/(1-z)], 


i.e.,  y  is  the  complex  natural  logarithm  of  the  integrand,  <  is  the  curvature  of 
the  osculatory  parabola  as  in  (9.8),  and  Xq  solves 

r(x0)  =  0.  (13.23) 

Subroutine  GETPAR  gets  the  parameters  of  integration  listed  above.  It  solves 
(13.23)  by  using  several  iterations  of  the  golden-section  algorithm  [2],  and  then 
it  switches  to  a  second-order  Newton's  method.  The  behavior  of  the  integration 
procedure  was  found  to  depend  somewhat  sensitively  on  the  accuracy  of  x0>  thus 
necessitating  the  extra  elaboration. 

Once  Xq  is  found,  the  second  derivative  <P"(Xq)  and  curvature  k  are 
determined  by  taking  differences.  The  curvature  k  is  found  directly  from  a 
difference  formula  assuming  local  second-order  behavior  (see  [5],  eq.  57).  The 
differences  in  the  independent  variable  are  taken  along  the  imaginary  direction 
in  order  to  avoid  complications  arising  from  the  proximity  of  the  point  (1 .0, 
0.0),  the  first  eigenvalue,  and  any  nearby  eigenvalues.  One  has  the  capability 
in  GETPAR  of  successively  refining  Ay  and  making  new  predictions  of  <  and  f " ( Xq ) 
by  extrapolation  to  Ay  =  0.  However,  the  zeroth  iteration  has  proved  adequate 
thus  far. 

(iii)  The  integration  procedure 

The  parabola  approximating  the  path  of  steepest  descent  is  parameterized  by 


using  y,  the  distance  from  the  real  axis,  as  the  parameter.  Thus,  one  has 


f  \  1  2  . 

z(y)  =  x0  2  Ky  +  iy’ 

since 

,  ,  12 
x(y)  =  x0  +  2  Ky 

on  the  osculatory  parabola.  One  capitalizes  on  the  fact  that  because  the  values 
of  the  integrand  above  Im  z  =  0  are  the  complex  conjugates  of  those  below,  one 
needs  only  to  follow  the  contour  for  positive  y  as  in  (9.10).  The  integration  is 
by  the  trapezoidal  rule,  which  is  rewritten  as  a  quadrature  formula  in  terms  of 
yk  =  kAy,  and  as  in  (9.11), 

Ay  =  O"(z0)]  /z. 

Finally,  one  has 

•  ;  I?  flzoH'  - 1  lss)4y 
00 

+  ^  f(Zj)(  l-iicy^  ) Ay } 

j=1 

with  f(z)  =  exp  [ V ( z ) ] .  The  summation  is  stopped  when  the  absolute  value  of  the 
summand  falls  below  epsilon  times  the  accumulated  sum,  where  epsilon  is  a  user- 
input  variable  of  the  order  of  10  .  The  integral  is  refined  by  repeating  it  for 
successively  halved  values  of  Ay.  Three  values  of  Qk  +  and  Qk~  are  printed  to 
permit  confidence  in  the  convergence  of  the  integral. 

(iv)  Validity  check 

The  algorithm  of  B6dard  [35]  was  used  to  compute  the  probabilities  pk  and 
thus  Qk*  by  summation  as  in  (9.1),  and  Qk+  =  1-Qk“,  for  a  single  Lorentz 
spectrum.  The  FRED  algorithm  was  reduced  nearly  to  a  single  spectrum  by  choosing 
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a  very  small  power  for  the  second  peak  and  by  displacing  it  only  slightly  from 
the  first.  Table  13.1  presents  the  results,  which  agree  within  the  ability  of 
FRED  to  simulate  a  single  spectrum.  Figure  13.1  is  a  plot  of  Table  13.1,  also 
showing  the  results  of  the  saddlepoint  approximation  (see  [5],  eqs.  (48)  and 
(49)),  which  is  not  guaranteed  to  be  accurate  near  the  mean.  The  results  of  the 
Fredholm-determinant  method  fall  so  close  to  those  of  Bedard  that  it  is  not 
possible  to  separate  them  on  a  plot. 

One  can  image  a  spectrum  consisting  of  two  displaced  Lorentz  spectra  whose 
width  parameters  are  permitted  to  approach  zero.  The  spectra  thus  become  delta 
functions,  and  this  permits  one  to  solve  the  integral  equation  for  D(z)  easily, 
and  thus  to  obtain  explicit  expressions  for  Qk  +  and  ,Qk~.  If  one  imposes  the 
condition  f^T  *  f2T,  and  insists  that  y^T  =  2-nk^  ,  y2 T  =  2irk2>  k1  and  k2  being 
distinct  integers,  then  the  formulas  are 


l=k 


Q  =1-0 
k  wk 


3  ”  noi  =  nfiT 


03.24) 


b  =  nQ2  =  nf2 T. 


A  comparison  of  results  from  this  formula  and  from  an  appropriately  chosen 
spectrum  in  program  FRED  is  given  in  Table  13.2.  Here  the  agreement  is  not  so 
precise  as  in  the  previous  example,  since  too  small  a  value  of  y^  =  y2  would 
cause  problems  in  the  algorithm.  However,  it  is  clear  that  there  is  no  serious 
disagreement . 

(v)  New  results 


In  this  section,  we  let  8  =  y^ '  be  the  frequency  offset  of  the  second 
Lorentz  spectrum  relative  to  the  first,  and  y1  and  y2  be  the  width  parameters  of 


k' 

k' 


Table  13.1 


Counting  statistic  Qk  +  for  a  single  Lorentz  spectrum, 
p-10,  nQ-5,  T-1.0,  calculated  by  the  BSdard  algorithm, 
the  Fredholm  determinant  algorithm,  and  the  saddlepoint  approximation. 


Q  +  Bedard 


Qi<+ ,  Fred 


Qk+,  Saddlepoint  approx’n 


2 

0.92939535 

0.92939536 

0.93404541 

4 

0.68036006 

0.68036009 

0.71836016 

6 

0.37859723 

0.37859724 

0.33516088 

8 

0.16831838 

0. 16831839 

0.151  29038 

10 

6.3351 307E-2 

6 . 3351 306E-2 

5.751 591 2E-2 

12 

2.1 133291E-2 

2. 1 133289E-2 

1 . 931 5245E-2 

14 

6.4594820E-3 

6.459*181  IE-3 

5.9339735E-3 

16 

1 .8521 763E-3 

1 .8521759E-3 

1 .7090362E-3 

18 

5.0661 61 5E-4 

5.0661 598E-4 

4 . 6941 9 1 3E-4 

20 

1  .3377483E-4 

1  .3377  477E-4 

1  .2446659E-4 
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Figure  13-1  Plot  of  Qk+  vs  counts,  k,  as  computed  by  the  Bedard  algorithm,  the 

Fredholm  determinant  method^  and  the  saddlepoint  approximation. 
Here,  p  =  10,  T  =  1.0,  and  n^  =  5.  For  the  Fredholm  determinant 

method,  p,  =  p0  =  10,  f,/f0  =  999. 


Table  13-2 


Counting  statistic  Qk  +  for  two  discrete  Lorentz  spectra, 
approximated  in  FRED  by  '  =  0.01,  y^"  =  10, 


nQ  =  5,  T  =  1.0,  a  =  3-75,  b 

=  1.25. 

k 

Qk+,  eq.  ( 13-24) 

Qk\  FRED 

2 

0.7805820 

0.77919718 

4 

0.5350658 

0.53318556 

6 

0 .3^8475^ 

0 .3147399 1 9 

8 

0.2218190 

0.221661126 

0 

0.1 396803 

0.14011555 

2 

8.7149901  E-2 

8.8187526E-2 

5.467138E -2 

5.5393858E-2 

6 

3.H11693E-2 

3.4762644E-2 

8 

2.1 27697E-2 

2.1805994E-2 

!0 


1.326526E-2 


1.36757H9E-2 
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the  spectra.  Figure  13-2  shows  the  plot  of  Qk  +  versus  k  for  two  spectra  with 
Pi  =  p2  =  10*  2  =  50,  and  values  of  f2/f.|  equal  to  1 ,  5,  and  10.  When  f1  =  f2, 
the  spectrum  is  maximally  bimodal,  and  as  f2/f.|  increases  it  becomes  unimodal,  so 
that  the  curve  approaches  that  of  Figure  13-1.  The  curves  diverge  perceptibly 
for  large  values  of  k.  Figure  13-3  is  the  same  type  of  plot,  but  here  we  let 
W1  =  v2  =  1  *0»  8  =  5.0,  T  =  1  .0,  and  f2/f.|  =1.5,  and  10.  This  is  a  case  where 
u-,T  =  p2T  is  small.  The  curv_o  have  the  same  general  behavior.  Since  the 
vertical  axis  is  drawn  on  a  probability  scale,  a  Gaussian  random  variable  would 
plot  as  a  straight  line,  and  obviously  the  variable  Qk  +  is  not  Gaussian. 

Actually,  these  plots  should  be  staircase  functions  continuous  from  the  right  at 
the  integers,  but  they  are  represented  here  as  curved  lines  for  simplicity. 

Figures  13-1*  and  13.5  are  cases  where  everything  about  the  two  Lorentz 
spectra  is  fixed  except  their  displacement  3  in  radian  frequency.  We  let  &/p1  = 
.001,  1.0,  and  10.0  on  each  figure.  For  Figure  13-2*,  =  u2  =  10»  fi  =  f 2  = 

0.5,  T  =  1.0,  and  nQ  =  5.  For  Figure  13-5,  we  have  =  p2  =  1  - 0 ,  and  otherwise 
the  variables  are  the  same.  In  these  cases,  the  increasingly  bimodal  character 
increases  Qk  +  for  a  given  k  when  k  >  nQ . 

The  coding  in  program  FEED  is  general  in  the  sense  that  it  does  not  care  how 
many  Lorentz  spectra  are  combined  to  make  up  the  power  spectral  density.  Only 
subroutine  FILPAR  and  some  dimension  statements  need  to  be  changed  to  handle  an 
arbitrary  number  of  spectra.  Thus  the  matrix  approach  to  finding  the  Fredholm 
determinant  and  the  counting  statistics  of  incoherent  narrowband  light  can  be 
applied  with  program  FRED  to  a  nearly  arbitrary  power  spectral  density. 

(c)  Coherent  plus  Chaotic  Light 

When  the  incident  light  consists  of  a  finite  number  of  sinusoidal  components 
in  the  presence  of  Gaussian  background  light  with  a  rational  spectral  density, 
the  moment -generating  function  h(z)  in  (13-7)  can  be  expressed  in  terms  of 

94 


0.0  3.0  6.C  9.' "  12.0  15.0  18.0 

COUNTS 


Figure  13.2  Plot  of  Qk+  vs.  counts,  k,  for  different  relative  powers  of  two 

component  Lorentz  spectra,  computed  by  the  Fredholm  determinant 
method.  Here,  y .  =  y „  =  10,  6  =  50,  T  -  1.0,  nn  =  5,  and  f,/f, 
5.  and  10.  1  2  0  21 


C001 


Figure  13.3  Plot  of  Qk+  vs.  counts,  k,  for  different  relative  powers  of  two 
component  Lorentz  spectra,  computed  by  the  Fredholm  determinant 
method.  Here,  y.  =  u2  =  1 ,  6  =  5 ,  T  =  1 .0,  n^  =  5 ,  and  f i  =  1 
5 ,  and  10.  '  ' 
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Figure  13.^  Plot  of  Qk+  vs.  counts,  k,  for  different  values  of  frequency  of 
separation  of  two  component  Lorentz  spectra,  computed  by  the 
Fredholm  determinant  method.  Here,  p.  =  p^  =  10»  =  =  0.5 

T  =  1.0,  nn  ■=  5,  and  B/pi  -  0.001,  1.6,  and  10.0. 


COUNTS 


Figure  13-5  Plot  of  Qk+  vs.  counts,  k,  for  different  values  of  frequency  of 

separation  of  two  component  Lorentz  spectra,  computed  by  the 
Fredholm  determinant  method.  Here,  ^  =  1'0,  fi  =  f2  = 

T  -  1.0,  nn  =  5,  and  8/y.  =  0.001,  1.0,  and  10.0. 


determinants  by  an  extension  of  the  analysis  in  Section  6.  The  complex  envelope 
of  the  normalized  signal  is  taken  to  have  the  form 


o(t)  «  T_l/2 


2Jctr  exp  (iv^t),  0  <  t  <  T, 


(13.25) 


in  which  the  vr  are  the  angular  frequencies  of  the  coherent  components  of  the 
light  and  the  ar  are  their  complex  amplitudes.  Hie  normalized  spectral  density 
?U)  of  the  background  light  can  be  written  in  the  form 


r(w)  =  2 


k-<  (U-U|<»)2  *  pk'2 


as  in  (6.5).  Its  normalized  autocovariance  function  is  then 


•y(t-u)  =  ^  ^  exp[-pk*(t-u)],  t  >  u, 

k  =  l 
n 

=  4>k*  expLuk(t-u)  ] ,  t  <  u, 


(13.26) 


k  =  1 


with  4>k  =  fk  +  igk,  yk  =  uk'  +  iyk".  The  spectral  density  r(w)  is  assumed  to 
have  only  simple  poles. 

We  first  obtain  an  expression  for  the  square  bracket  multiplying  nQ(z-1)  in 
(13-7);  we  denote  it  by  J.  In  terms  of  the  eigenfunctions  and  eigenvalues  of  the 
integral  equation  (13*9),  it  is  [25] 


J(0  = 


X""*  l°iJ  f 

2^  i  +  IT'  °k  =  J  fk*(t)  o(t)  dt’ 

K  J^ 


o  k |  =  1i  £  =  no( 1 _z) ’ 


(13-27) 
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Now  define  the  signal  q(t;£)  as  the  solution  of  the  integral  equation 


(t)  =  q(t;0  +  C T  /  Y(t-u)  q(u;£)  du,  0  <  t  <  T.  (13.28) 


It  is  not  difficult  to  show  from  the  orthonormality  of  the  eigenfunctions  fk(t) 
of  (13.9)  that 


J(0  =  /  o* ( t )  q(t;£)  dt. 


(13.29) 


The  integral  equation  (13*28)  has  a  solution  of  the  form 


q(t;£)  =  q00(t;U  dk  exp  (BRt), 


(13-30) 


in  which  as  in  Sec.  6  the  Bk's  are  the  roots  of  (13.19),  and 


-y  v-^  r 


a  exp(iv  t) 


*(v)  =  1  +  £T  f(v).  (13-31 ) 


Substituting  (13. 30)  into  (13.28),  one  finds  that  the  coefficients  dk  have  the 


dK  ‘  T  7/Ar 


(13-32) 


with  the  dkr  the  solutions  of  the  2n  linear  simultaneous  equations 


d,  exp(  6,  T)  exp(  iv  T) 
kr  k  r 


8,  -  u 
k 


<Kv„)(pm-iv„)  ’ 
r  m  r 


1  <  m  <  n, 


8^  -  u  ♦  (  v  )  ( it  -iv)  ’ 
k  m  r  m  r 


n+1  <  m  <  2n,  (13*33) 
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'Vs 


(13.38) 


%  %  V 


~sk 


exp[(6R-ivs)T] 
""  (ek  '  ivs)Y 


With  (13*35)  this  becomes 


--EE 


2n  2n 


2n  2n 


J2  - 


EkgkmFm’ 


r  s 


k  =  1  m  =  1 


k  =  1  m  =  1 


o*e  ,  ,  F  =  >  ct  f  /*(v  ). 
3  sk  m  /  j  r  mr  r 


(13.39) 


By  defining  the  (n+1)  x  (n  +  1)  augmented  matrix  G  as 

'0  E^ 


G  = 


F  G 


(13.^0) 


where  E  is  the  2n-element  row  vector  of  the  Ek's  and  F  is  the  2n-element  column 


vector  of  the  Fm's,  we  can  write 


J2  =  -  det  G. 


03. **1 ) 


The  Fredholm  determinant  D(z)  in  (13.7)  is 

2n  n 


D(z)  =  exp  SkT  ^  um*T) 


det  G(T) 
det  G(  0 )  ’ 


(1  3.-42) 


k  =  1 


m  =  1 


where  G(0)  is  the  2n  x  2n  matrix  obtained  by  setting  T  =  0  in  the  elements  of  the 
matrix  G  =  G(T).  The  probability-generating  function  of  the  photoelectron  counts 
is,  as  in  (13-7), 


h(z)  =  [D(z ) ]  1  exptns(z-1 )](J1  +  J2 ) ] , 


(13.^3) 


E  =  nQ(1 -z)  ; 

Jl  and  J2  depend  on  5  through  the  Bk's  and  4>(v)  defined  in  (13-31)  and  (13.19); 
4>(-igk)  =0,  1  <  k  <  2n. 

For  a  single  sinusoidal  component,  0^  =  1 ,  v1  =  v,  J1  =  T/fl>(v), 
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-  ’  *-  %  A 


.-a  i.jr*...*; 


F,  =  - ,  n  +  1  <  k  <  2n ,  y,  =  y*  .  (13.44) 

k  y.  -  iv  -  kk  pk-n 

k 

The  expression  in  (13. 11*)  is  a  special  case  of  (13.43)  and  (13*^^)  for  n  =  1. 

This  formulation  enables  the  probability  generating  function  h(z)  to  be 
evaluated  entirely  by  finite  operations  for  all  values  of  z  involved  both  in 
locating  the  saddlepoints  of  the  integrands  of  (9.4)  and  (9.5)  and  in  evaluating 
those  integrals  for  the  cumulative  probabilities  by  numerical  quadrature  in  the 
complex  plane.  The  evaluation  entails  solving  the  algebraic  equation  of  degree 
2n  with  complex  coefficients  arising  from  (13.19)  and  calculating  the 
determinants  in  (13*41)  and  (13-42);  standard  computer  routines  can  be 
employed.  Any  finite  number  of  sinusoidal  components  can  be  handled. 

If  the  coherent  signal  S(t)  is  a  pulse  confined  to  the  counting  interval 
(0,  T),  the  complex  amplitudes  ap  in  (13.25)  can  be  determined  by  numerical 
Fourier  transformation  of  the  normalized  signal  o(t);  a  large  enough  number  of 
Fourier  components  must  be  utilized  to  provide  adequate  accuracy  in  the 
cumulative  distribution.  The  term  J1  in  (13-37)  and  (13.43)  becomes 


(13.45) 


with  $(•)  as  in  (13-31)  and  with  Km)  the  Fourier  transform  of  the  normalized 
signal  o(t);  this  is  just  the  first  integral  in  the  Toeplitz  approximation  in 
(13*20).  The  terms  Ek  and  Fm  in  (13-39)  become 
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F 

m 


expCB^t)  dt/T  = 

[eluT]_ 


$(u))(p  -ioj) 
m 


dui 
2  TT  ’ 


T  (Z(-i6k*) )*. 


03. 46) 


in  which 

r  iwT-,  iwT  . 

Le  ]  =  e  ,  1  <  m  <  n, 

m  - 

=  1 ,  n+1  <  m  <  2n. 

The  details  of  applying  these  new  methods  to  computing  the  photoelectron  counting 


distribution  remain  to  be  elaborated. 


14.  Significance  Probabilities  for  Rank  Tests 
The  method  of  saddlepoint  integration  described  in  Sec.  9  has  been  applied 
to  computing  the  distributions  of  the  Wilcoxon  signed-rank  and  ranksum  statistics 
[49]  and  the  Jonckheere  statistic  [50].  Tables  of  the  distributions  of  the 
Wilcoxon  statistics  and  of  their  percentage  points  take  up  a  great  many  pages 
because  of  the  numbers  of  parameters  involved  [49],  and  the  distribution  of  the 
Jonckheere  statistic,  which  is  a  generalization  of  the  ranksum  statistic  and 
depends  on  even  more  parameters,  seems  to  have  been  tabulated  only  for  a  few  ve^y 
small  sample  sizes  [50],  Although  recurrent  relations  exist  whereby  these 
distr ibutions  can  be  computed,  the  algorithms  are  extremely  time-consuming  and 
require  extensive  computer  memory  when  the  sample  sizes  are  large.  By 
saddlepoint  integration  one  can  compute  the  significance  probabilities  directly 
and  accurately  for  sample  sizes  beyond  those  at  which  the  recurrent  algorithms 
necome  mfeasiDie.  Tne  saddlepoint-integration  method,  furthermore,  requires 
s'.  :.".ng  only  a  dozen  or  so  numbers  utilized  in  the  computation  ,  in  contrast  to 
t-.-  r.  u.d-  eds  cr  thojsands  that  may  need  to  be  stored  in  the  recurrent 

1  3: gn-c  :-Ran-:  Statistic 

7r-  one-sample  Wilcoxon  test,  or  "signed-rank  test,"  decides  between  a  null 
nyp  .  tres :s  H„  tnat  M  independent  data  x  ,  x  ,  ...,  x  possess  a  common 

J  1  c.  M 

dist” ibuti  ?n  symmetrical  about  0  and  an  alternative  hypothesis  H1  that  the  data 

a-fe  stochastically  larger.  The  data  are  first  ordered  according  to  their 

absolute  values.  If  r.  denotes  the  rank  of  x^  in  that  ordering,  one  forms  the 

statistic 

M 

2Zri  “'V' 

i  =  1 


W 


(14.1  ) 


where  U(*)  is  the  unit  step-function.  When  a  given  set  of  Jcta  yields  a  value  W 


for  this  statistic,  one  quotes  the  significance  probability 

P  =  Pr  (W  >  W^|  Hq) 


that  a  value  of  W  as  large  as  W^  or  larger  could  occur  under  the  null  hypothesis 
and  one  states  that  the  null  hypothesis  is  rejected  at  probability  level  P. 
Detailed  explanations  of  this  test  and  its  two-sided  extension  can  be  found  in 
texts  on  statistics. 

Because  the  distribution  of  W  under  hypothesis  HQ  is  symmetric  about  its 


mean 


W  =  E (W [  Hq)  =  M(M+1  )/4,  (1  4.2) 

this  probability  is  the  same  as  the  probability 

P  =  Pr  (W  <  WQ|  Hq),  WQ  -1  M(M+1 )  -  W£. 

We  need,  therefore,  the  cumulative  distribution  of  W,  which  by  (9.4)  is  given  oy 
the  contour  integral 


q^  =  Pr  ( W  <  k  )  = 


(z)  dz 
2iri 


(14.3) 


in  terms  of  the  p.g.f.  h(z)  of  the  signed-rank  statistic  [49] 


h(z) 


M 

IW 


(14.4) 


j=1 


Here  C  is  initially  a  circle  concentric  with  the  origin  and  of  radius  less  than 

1  . 


It  has  been  found  simplest  to  deform  the  contour  C  into  a  vertical  straight 


line  through  the  saddlepoint  Xq  of  the  integrand  of  (14.3)  lying  on  the  Re  z- 
axis,  0  <  Xq  <  1,  completed  by  the  portion  of  the  unit  circle  lying  to  the  left 
of  that  straight  line,  whereupon 


q„  =  Re 


/ 


-(k+1 )  . 

z  h(z) 


dy 

IT 


^  +  I 


2  =  X  , 


iy. 


2  V 

■  <'-«o  >  • 


(14.5) 


The  second  term 


I  = 


/•IT  M 

/n 

en  j=i 


cos 


(j  J0) 


Re  [ 


exp[i(W-k)6]i  de 


l-exp(ie) 


1 

J  7T  ’ 


e„  = 


-i 

cos  x 


o’ 


(14.6) 


represents  the  integral  along  the  unit  circle,  on  which  z  =  exp(ie). 

For  the  large  values  of  M  for  which  this  method  is  intended,  the  term  Ic  is 
neglibible,  and  sufficient  accuracy  is  achieved  by  integrating  only  the  first 
term  in  (14.5),  The  product  of  cosines  in  (14.6)  possesses  M(M+1)/4  zeros,  some 
multiple,  in  0  <  e  <  it.  Between  any  pair  of  these  zeros,  which  lie  close 
together  when  M  is  large,  the  oscillatory  integrand  perforce  remains  very  small, 
and  the  value  of  Ic  itself  must  be  very  small  in  comparison  with  the  first  term 
in  (14.5).  Our  computations  have  borne  out  the  insignificance  of  the  term  Ic 
when  M  is  greater  than  about  20,  as  manifested  by  the  relative  errors  to  be 
exhibited  in  Table  14.1. 

The  saddlepoint  xQ  is  the  root  of  the  equation 


where 


y  *  (x ) 


M  .  , 

J=1 


k+1 

x 


(14.7) 


M 

f (z)  °  ^  ]  In  (--*-  z  )  -  (k  +  1)  In  z  -  In  (1-z) 

J-1 


1C'7 


(14.8) 


d  .  fiiTi.W.i  ■'  Wl  A.  A.  W- A.  .f- *-» Ik .tv A' Al ,*>  -*»  1. -*»VW‘«Y .  ^  ^  » 
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Table  14.1 


Signed  Rank  Statistic 


k 

% 

Number 
of  steps 

Relative 

error 

21 

4 .24988 ( -4 ) 

15 

-8.2(-4) 

4.24993(~4) 

30 

-8.1 (-4) 

4 . 24994 (-4 ) 

60 

-8.K-4) 

43 

9 .61 708 (-3) 

11 

2. 4 (-5) 

9.6l7l0(-3) 

22 

2 .6  ( -5 ) 

9 .6171 4 (-3) 

43 

3 . 0  ( —  5 ) 

81 

0.194199 

5 

5.7(-5) 

0.194188 

10 

-2. 2(-6 ) 

0.194187 

20 

-4. 3 (-6) 

78 

4.75886(-4) 

11 

-1 .6 (-6) 

4.75889(-4) 

22 

3 . 3 ( -6 ) 

4.75890(-4) 

44 

5.7 (-5) 

120 

9.83041 (-3) 

6 

6 . 0 ( -7 ) 

9.83043(-3) 

12 

2. 5 (-6 ) 

9 .83043 (~3 ) 

24 

2 .5 ( -6 ) 

190 

0.196536 

6 

5.5 (-5) 

0.196525 

1  1 

-6 . 4 (-8) 

0.196525 

21 

-5.3C-7) 

172 

4 .968 1 7 ( -4 ) 

7 

— 8 . 5  ( -  7 ) 

4 . 96818 (-3) 

14 

8.8(-7) 

4.9  68 1 8 ( - 3 ) 

28 

1  .7 ( -6 ) 

238 

9 . 9939  7 ( ~  3 ) 

6 

5.6 (-7) 

9 .99396 (-3 ) 

1 1 

-6 .7 (-7 ) 

9 .99395 (-3) 

22 

-1 .5 (-6) 

346 

0.198723 

6 

5  .5  ( ~5 ) 

0.198712 

1 1 

-7.6 (-8) 

0.198712 

21 

2 .5 (-7 ) 

t  V 


k 

<»k 

Number 

of  steps 

Relative 

error 

689 

4.96333(-10 

6 

3 .7  C  —8 ) 

4.96333(-10 

12 

-9.4(-8) 

4.96333(-4) 

23 

-2.7C-7) 

846 

9.91121 (-3) 

6 

7 .8 (-7 ) 

9.91120C-3) 

1 1 

-2.6(-8) 

9.91 1 20 (-3) 

22 

3. 6 ( —8) 

1097 

0.199230 

6 

5.5(-5) 

0.199219 

11 

1 . 4  ( —  7 ) 

0.199219 

22 

3.6 (-7 ) 

1578 

4.94497(-4) 

6 

4.7 (-8) 

4.94497C-4) 

1 1 

8 . 9  ( -8 ) 

4.94497C-4) 

22 

2. 9 (-7) 

1850 

9 . 965  4  8 ( —  3 ) 

6 

9.1  (-7) 

9 .96547 (-3) 

1 1 

3.5 (-8) 

9 .96547 (-3) 

22 

1  . 0  ( —  7 ) 

2278 

0.199166 

6 

5.6 (-5 ) 

0.199155 

12 

1 .4 ( -8 ) 

0.199155 

23 

1  -  4 (-7 ) 

show  how  the  integrals  evaluated  by  the  trapezoidal  rule  in  double-precision 
arithmetic  converge  as  the  step  size  is  halved.  We  took  e  =  10“3  in  truncating 
the  numerical  integration.  The  relative  error  tabulated  there  is  the  difference 
between  the  probability  as  calculated  by  saddlepoint  integration  and  the  exact 
probability,  divided  by  the  latter.  The  exact  probability  was  computed  by  the 
recurrent  method  of  Mann  and  Whitney  [51]  carried  out  in  double-precision 
arithmetic.  These  computations  were  executed  on  a  VAX  11/780. 

In  order  to  compare  the  time  required  by  the  saddlepoint  integration  method 
with  that  required  by  the  recurrent  algorithm  [51]  we  programmed  both  in  Basic  on 
an  HP-71B  handheld  computer.  Results  are  exhibited  in  Table  14.2.  The  recurrent 
algorithm  may  need  to  compute  and  store  as  many  as  W  =  M(M+1)/4  values  of  the 
probabilities  pk  =  Pr  (W  =  k),  and  the  time  TR  listed  in  Table  14.2  is  that 
needed  to  do  so.  The  time  to  sum  those  probabilities  to  obtain  the  cumulative 
probability  qk  is  not  included.  The  times  in  the  column  labeled  Tj  are  those  for 
the  entire  computation  of  qk  by  saddlepoint  integration.  Again  we  took  c  =  10~^; 
Nst  is  the  number  of  steps  in  the  numerical  quadrature. 

Once  the  number  M  of  samples  exceeds  thirty  or  so,  the  saddlepoint 
integration  method  is  decidedly  faster  than  the  recurrent  algorithm.  It  should 
be  kept  in  mind  that  the  statistician  usually  needs  only  a  single  significance 
probability  Pr  (W  <  WQ)  in  assessing  a  set  of  experimental  data.  The  recurrent 
algorithm  must  compute  and  store  the  probabilities  pk  for  all  values  of  k  up  to 
and  including  k  =  WQ  and  sum  them,  and  if  Wq  is  near  the  mean  W,  this  may  take  a 
long  time  and  require  storing  many  numbers. 

(b)  The  Ranksum  Statistic 

In  the  two-sample  Wilcoxon  test,  or  "rank -sum  test,"  one  arranges  the  N+M 
data  x1 ,  x^,  ....  xM,  y1 ,  y^,  ....  y^  in  ascending  order.  Let  r  denote  the  rank 
of  the  datum  x^  in  that  arrangement.  One  forms  the  rank-sum  statistic 
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Table  14.2 

Signed-rank  Statistic 


k  Nst  TjCaec.)  Qk 


Relative 

error 


M 

=  25, 

T  =  67.6  sec. 

n 

45 

14 

78.0 

4.55944(-4) 

3-11 (~6) 

76 

8 

41  .4 

9.37229(-3) 

6. 39 (-6) 

130 

6 

29.8 

1 . 97867 (-1 ) 

— 5 . 68 ( — 5 ) 

M 

=  30, 

Td  =  117.6  sec. 

n 

78 

12 

77.0 

4.75860(-4) 

5.81 (-5) 

120 

7 

43.2 

9 .8301 1 (-3) 

3.02(-5) 

190 

7 

39.4 

1 .96541 (-1 ) 

-8.21 (-5) 

M 

=  40, 

Td  =  281 .6  sec. 

K 

172 

8 

70.9 

4.96793(-4) 

4.93C— 5) 

238 

7 

53.4 

9 .99337 (-3) 

5 . 99 ( —5 ) 

346 

7 

50.8 

1 .98727 (-1 ) 

-7.57 (-5) 

M 

=  50, 

T  =  554.0  sec. 

R 

304 

8 

83.9 

4.93541 (-4) 

4.06C-5) 

397 

7 

66  .0 

9.79352(-3) 

3  - 76 ( -5 ) 

549 

7 

63-3 

1 .99441 (-1 ) 

-6.82C-5) 

If  for  a  particular  set  of  data  this  yields  the  value  Vq,  one  quotes  the 
significance  probability  P  =  Pr  (V  >  Vq)Hq)  that  as  large  a  value  of  V  as  VQ,  or 
a  larger  value,  could  arise  under  the  null  hypothesis  HQ  that  the  x's  have  the 
same  distribution  as  the  y's;  and  one  reports  that  the  null  hypothesis  is 
rejected  at  probability  level  P. 

The  maximum  possible  value  of  V  is  NM  +  ^  M(M+1).  The  significance 
probability  P  is  given  by 

P  =  Pr  (V  >  VQ|  HQ)  =  Pr  (U  <  UJ  HQ)  ,  (14.14) 

where 

U  =  NM  +  ~  M(M  +  1 )  -  V,  UQ  =  NM  +  |  M(M+1 )  -  VQ .  (14.15) 

The  random  variable  U  has  the  following  significance.  Consider  all  ,,  ..C., 

N+M  M 

arrangements  of  M  x's  and  N  y's;  under  HQ  these  are  equally  li<ely.  Then  in  a 
randomly  selected  one  of  these  arrangements  U  is  the  number  of  times  that  an  x 
precedes  a  y.  The  value  of  U  ranges  from  0  to  NM,  and  its  distribution  is 
symmetrical  about  its  mean 

U  =  E(Uj  Hq)  =  MN/2.  (14.16) 

We  treat  the  computation  of  the  cumulative  probabilities 

qk  =  Pr  (U  <  k|H0)  (14.17) 

for  0  <  k  <  MN/2,  from  which  the  significance  probability  P  =  Pr  (V  >  VQ|  HQ)  can 
be  calculated  by  (14.14)  and  (14.15).  For  the  sake  of  efficiency  one  will  always 
take  M  as  the  smaller  of  the  two  sample  sizes,  M  £  N,  interchanging  the  labels  x 


and  y  on  the  data  if  necessary. 


The  probability  generating  function  of  the  random  variable  U  is  [49] 


h(z) 


(n+mV 


-1 


(14.18) 


Again  we  can  write  the  probability  qk  as  the  contour  integral  (14.3),  and  again 
we  pick  a  closed  contour  C  surrounding  the  origin,  symmetrical  about  the  Re  z- 
axis,  excluding  the  point  z  =  1,  and  passing  vertically  through  the 
saddlepoint  z  =  xQ  of  the  integrand  lying  in  0  <  x  =  Re  z  <  1.  The  phase  of  the 
integrand  is  now 


M 

¥(z)  =  'y  ^[ln  (zN+j-1)  -  ln(z^-.l)] 
j=1 

+  In  c  -  (k  +  1)  In  z  -  In  (1-z), 


(14.19) 


and  the  saddlepoint  xQ  is  the  root  of  the  equation 

M 

r  (z)  = 

j  =  1 


(N+j )  z 


N+j  -1 


zN+j-1 


-  i£l]  -  ill .  _L  .  o 
ZJ-1  z  z~' 


(14.20) 


lying  in  that  interval.  As  before,  we  integrate  instead  along  the  straight 

vertical  line  segment  passing  through  the  saddlepoint  z  =  xQ  and  lying  within  the 

unit  circle,  and  we  integrate  thence  counterclockwise  around  the  unit  circle. 

The  contour  integral  (14.3)  takes  on  the  same  form  as  in  (14.5),  except  that  now 

the  integral  I  along  the  unit  circle  to  the  left  of  the  line  Re  z  =  xn  is 
c  0 

written 


.1  H  1 

/  n— 

J  7  sin 
J=1 


j  (N+j)e 


Re  [ 


exp 


2  Je 


i( U-k ) 6  I  d9 
exp(i9)  it  * 


-1 

cos  x 


(14.21) 


The  same  argument  as  in  part  (a)  leads  us  to  expect  that  when  M  and  N  are 


large,  the  term  I  can  be  neglected.  The  product  in  the  integrand  of  (14.21)  nas 
MN/2  zeros  —  some  of  them  multiple  —  lying  in  0  <  6  <  ir,  and  the  integrand, 
which  is  oscillatory,  perforce  remains  small  between  them.  The  principal 
contribution  to  comes  from  the  first  term  of  (14.5),  with  h(z)  as  in  (14.18), 
and  in  that  term  it  arises  from  the  neighborhood  of  the  saddlepoint  z  =  xQ.  The 
larger  M,  N,  and  k,  the  more  rapidly  the  integrand  in  that  term  drops  off  from 
its  peak  value  at  y  =  0. 

The  first  term  of  (14.5),  with  (14.18),  is  integrated  by  the  trapezoidal 

rule,  and  the  integration  is  stopped  when  the  absolute  value  of  the  integrand 

falls  below  e  times  the  accumulated  sum  times  the  step  size  Ay,  provided  y  <  Y; 

-3 

again  we  took  e  =  10  .  The  step  size  Ay  is  chosen  as  in  (14.13),  with  the 

second  derivative  4"'(x)  of  the  phase  again  approximated  as  in  (14.9),  where  now 

p  =  Var  U  -  U  =  MN(M+N~5)/12.  (14.22) 

The  saddlepoint  is  calculated  as  in  part  (a)  by  solving  (14.20)  by  Newton's 
method  (1.16),  the  initial  trial  value  being  determined  by  (14.12)  with  W 
replaced  by  U  =  MN/2  and  p  by  (14.22). 

In  Table  14.3  we  show  how  the  values  of  the  probability  qk  calculated  by 
numerical  integration  converge  as  the  step  size  is  halved.  Exact  values  of  the 
probabilities  were  computed  by  Harding's  algorithm  [52]  carried  out  in  double¬ 
precision  arithmetic,  and  the  relative  errors  were  determined  as  in  part  (a). 

The  numerical  quadrature  was  terminated  as  before,  with  e  =  10-^.  Good  accuracy 
is  attained  with  only  a  few  steps  of  numerical  integration. 

Table  14.4  compares  the  times  required  by  an  HP-71 B  computer  to  carry  out 
the  numerical -quadrature  algorithm  and  Harding's.  The  time  TR  is  that  occupied 
by  the  latter  when  computing  and  storing  all  the  probabilities  Pk>  0  <  k  <  U,  as 

?  r 

it  potentially  must  do;  TR  increases  as  M  N  as  M  and  N  increase  l52j.  The  times 


Table  14.3 


Ranksum  Statistic 


k 

% 

Number 
of  steps 

37 

4 . 983  38 ( — 4 ) 

7 

4.98342(-4) 

14 

4 . 98343  C  —4 ) 

28 

57 

9.27685(-3) 

6 

9.27684(-3) 

1 1 

9.27682(-3) 

21 

92 

0.194630 

5 

0.194617 

10 

0.194617 

19 

93 

4 . 88862 ( -  4 ) 

7 

4.88863(-4) 

13 

4 .88862 ( - 4 ) 

25 

1  29 

9 . 4564  2 ( - 3 ) 

6 

9.45641 (-3) 

1 1 

9 . 456  40 ( -3 ) 

21 

190 

0.198499 

5 

0.198487 

10 

0.198487 

20 

1  38 

4.75283C-4 ) 

6 

4.75283(-4) 

12 

4.75282(-4) 

23 

183 

9 . 48660 ( - 3 ) 

5 

9 .48659 ( — 3 ) 

1  0 

9.48660(-3) 

20 

257 

0.196244 

6 

0.196232 

1 1 

0.196232 

21 

Relative 

error 

-2 .6 (-5 ) 
-1 . 9 (-5) 
-1 .6 (-5) 

-1 . 8 ( -6 ) 
-2.K-6) 
-5.1 (-6) 

6 . 6 ( —5 ) 
1.81-7) 
2 . 4 ( —6  ) 

— 8 . 4  ( — 7 ) 
2.0(-6) 
-4.1 (-7) 

6.8C-7) 
-2 . 3 ( -6 ) 
-1  .7 (-6) 

6.1 (-5) 
4 .6 ( -7  ) 
1 . 4 (-6) 

3 . 9 ( —7 ) 
-7.2(-8) 
-2. 1 (-6) 

6 . 3  C  -7 ) 
—2 . 3  C  —7  ) 
6 . 0  ( —7 ) 

6.0 (-5) 
2.0(-8) 
4 .9  ( -7  ) 
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Table  14 

.4 

Ranksum  Statistic 

Tj  (sec . ) 

10,  N  =  40,  Tr 

\ 

=  120.2  sec. 

Relative 

error 

52.8 

4 . 59224 ( -4 ) 

2 . 4  4  (  -  6  } 

30.3 

9 . 4 1 982 ( -3 ) 

-3.50' -t 

25.4 

1  .97525 (-1  ) 

5.63;-5, 

20,  N  =  20,  Tr 

=  247.5  sec. 

50.9 

4 .6671 0 ( -4 ) 

2. 29,-6, 

51.1 

9 .76839 (-3) 

5.521-7) 

44.9 

1  . 99 1 6  4 ( -1  ) 

6.37  : -5, 

15,  N  =  45,  Tr 

=  316.7  sec. 

55.1 

4  .67736  ( -4 ) 

-2. 94, -7, 

42.6 

9 .551 36 ( -3 ) 

6.281-7; 

40.0 

1  . 98  6  6 2  (  - 1  ) 

5.751-5/ 

20,  N  =  40,  Td 

n 

=  510.4  sec. 

57.5 

4 . 804 1 0 (-4 ) 

2.291-7) 

54.3 

9.90422(-3) 

8.531-7) 

51  .0 

1  .97985 ( -1  ) 

5.861-5) 

30,  N  =  30  Tr 

=  870.1  sec. 

72.7 

4 . 87  380 ( -4 ) 

6.321-8) 

68.6 

9.95029 (-3) 

9.221-7) 

72.8 

1 . 9901 4 ( -1 ) 

5.991-5) 

Tj  for  computing  the  probabilities  qk  by  saddlepoint  integration  are  considerably 
shorter,  and  they  are  proportional  roughly  to  the  smaller  of  M  and  N.  It 
requires  storing  only  a  dozen  or  so  numbers;  the  recursive  Harding  algoritnm 
requires  storing  as  many  as  11  =  MN/2,  and  this  might  well  strain  the  capacity  of 
a  small  computer. 

(  c )  The  Joncwheere  Statistic 

The  lonckheere  test  extends  the  rank-sum  test  to  an  arbitrary  number  s  >  2 
of  trie  samples  of  sizes  M1  ,  M  [50].  One  has  s  groups  of  independent 

data 


(Xu  >  x ( 1  'J  x(i))  i  =  i  2  s 

l 

Wit.n  U.  ,  denoting  trie  number  of  pairs  (a,  b)  for  which  x(-i  ^  <  x^  \  i  <  j  ,  the 
1  d  D  7 


;est  statistic  is 


s  j-1 

£  E.v 

j  "2  i=1 


(14.23; 


For  s  =  2,  M  =  M,  M?  =  U  is  tne  ranksum  statistic  treated  in  part  (b).  We 
are  concerned  with  tne  distribution  of  U  under  the  null  hypothesis  that  all  the 
data  xw'  have  a  common  probability  distribution. 

The  probability  generating  function  of  U  under  the  null  hypothesis  H  is 


5 

c  =  j~J  NM/t^S .  (14.24) 

m=1 


with 


If  =  M,  +  M-  +  ...  +  M 
t  1  2  S 


0  4.25) 


the  total  number  of  data.  Harding  [52]  has  shown  how  the  probabilities 

=  Pr  (U  =  nj  H) 

can  be  computed  by  s  -  1  repetitions  of  his  minimal-storage  algorithm  for  the 
ranksum  statistic. 

Saddlepoint  integration  can  conveniently  be  applied  to  computing  the 
significance  probabilities 

k 

qR  =  ^  Pn  =  Pr  (U  <  k |  J|).  (14.26) 

n=0 

It  suffices  to  consider  the  range  0  <  k  <  U,  where 

s 

U  =  J  (Nt2  -  ^  M.2)  (14.27) 

J-1 

is  the  mean  of  the  distribution,  about  which  it  is  symmetrical, 

qk  =  1  q2U-k ' 

The  probability  is  again  written  as  the  contour  integral  04.3). 

For  the  sake  of  fast  computation  we  rewrite  h(z)  in  terms  of  the  numbers 
p.,  pof  . ..,  p  of  different  sample  sizes,  of  which  there  are  p.  Let  r.  be  the 

1  2  p  j 

multiplicity  of  size  p j  .  If,  for  example,  s  =  5  and 

M1  =  M2  =  5,  M3  =  =  M5  =  10, 

then  p  =  2,  p1  =  5,  r 1  =  2,  p2  =  1 0,  r2  =  3.  If  all  s  samples  have  the  same  size 
M,  then  p  =  1,  p1  =  M,  =  s.  Now 

1 20 


.a//  /  V .  * 


S' 

Nt  =  rjyj  ’ 


(14.28) 


the  mean  is 


I  (Nt2  -Erj“j2)’ 


(14.29) 


and  the  p.g.f.  is 


/  <  fn  v 

(1-z  ) 


p  Mj 


r . 
"M  J 


j=1  m  =  1 


(14.30) 


For  later  convenience  the  sample  sizes  y.  are  arranged  in  ascending  order, 

J 


h  <  w2  <  •••  <  V 


The  phase  of  the  integrand  in  (14.3)  is 


t(z)  =  In  c  +  )  ln(1-z  ) 


CO 


P  yJ 


r^  ^  ln(1-zm)  -  (k+1)ln  z 
j=1  m  =  1 


-  In ( 1  -z) .  (14.3D 


Its  derivative 


(z)  =  z 


_  «  I “2  «  «  _  *  1  2 


j=1  m-1  . 


is  again  set  to  zero,  and  the  resulting  equation  is  solved  for  the  saddlepoint  xQ 
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lying  in  (0,  1).  Newton's  method  (1.16)  is  used,  in  which  the  second  derivative 
^"(x)  is  approximated  as  in  (14.9),  where  now 


P 

5  ■  k  K3  -]£vj3 ' 3051  (’“-33) 

j=1 

The  integration  is  again  taken  over  a  straight  vertical  contour,  completed 
around  the  part  of  the  unit  circle  lying  to  the  left  of  it,  as  in  (14.5). 

Because  the  p.g.f.  h(z)  is  a  product  of  p.g.f.'s  of  the  type  in  (14.18),  h(z) 
will  have  a  large  number  of  zeros  on  the  unit  circle  |z|  =  1  when  N  >>  1,  and 
the  contribution  of  the  term  1  will  be  small.  The  integral  along  the  vertical 
path  is  approximated  by  the  trapezoidal  rule  with  a  step-size  Ay  given  by 
(14.13),  with  (14.9)  and  (14.33).  The  integration  is  stopped  when  the  absolute 
value  of  the  integrand  falls  below  e Ay  times  the  accumulated  trapezoidal  sum. 
Provided  this  occurs  before  the  point  z  =  xQ  +  iy  reaches  the  unit  circle,  the 

-3 

term  Ic  can  be  neglected.  We  found  e  =  10  sufficiently  small  for  accurate 
evaluation  of  q,  . 

The  bracketed  portion  of  (14.32)  is  easily  programmed.  One  accumulates  the 

sum 

Nt 

m(1-z  m)  1 

m  =  1 


in 


the  first  term,  picking  off  the  values  of  the  sums 


E 

m  =  1 


m(  1 


-tn .  -1 
-z  ) 


(14.34) 


for  the  second  term  when  m  passes  through  successive  values  of  p..  The  p.g.f. 

0 

h(z)  in  (14.30)  is  computed  in  a  similar  manner.  As  the  product 


/  4  m  \ 
(1-2  ) 


in  its  numerator  is  being  accumulated,  the  products 


, ,  m. 
(1-z  ) 


in  its  denominator  are  picked  off  when  m  passes  through  the  successive  sample 
sizes  \i.,  1  <  j  <  p.  Tne  time  needed  by  the  saddlepoint  integration  method  is 
roughly  proportional  to  N  .  Written  to  accommodate  a  variety  of  sample  sizes  and 
numbers  s  of  samples,  our  program  is  not  so  efficient  for  s  =  2  as  the  one  based 
on  the  equations  given  in  part  (b). 

Table  14.5  compares  the  performance  of  this  method  with  that  of  Harding's 
algorithm  [52].  The  basis  of  comparison  is  the  same  as  in  Table  14.4. 

Individual  significance  probabilities  qk  are  computed  in  markedly  less  than  the 
time  T  potentially  required  by  the  recursive  algorithm,  even  for  modest  sample 

n 

sizes.  Harding’s  method  must  allow  for  the  possibility  of  having  to  sto^e  Tu'' 
numbers,  with  U  given  in  (14.27);  the  integration  method  occupies  only  a  small, 
fixed  amount  of  computer  storage. 


Table  14.5 

Jonckheere  Statistic 

N  V 

St 

T  (sec . ) 

qk 

Relative 

1 

error 

a  -  3.  M1  - 

5,  M2  =  10, 

M3  =  15,  Tr  =  127.9  sec. 

6 

7 

60.5 

^  .8931 3C-1*) 

2.59 ( -7 ) 

8 

6 

48.6 

9  .55420 (-3) 

1 . 00 (-6 ) 

5 

6 

48.8 

1 .9601 1 (-1 ) 

6 .78  C —5 ) 

s  =  3,  M 

‘l  =  M2  M3 

=  10,  Td  =  189.1  sec. 

n 

5 

7 

56.6 

4.53680(-4) 

2.29(-7) 

8 

6 

45.6 

9 .079 1 3 ( -3 ) 

1.18C-6) 

7 

6 

45.5 

1 .98204(-1 ) 

6.95C-5) 

s  =  4,  M1 

=  m2  =  m3  = 

=  10,  Tr  =  591.2  sec. 

5 

6 

62.2 

4.56798(-4) 

2.42(-8) 

■3 

6 

58.9 

9 .52746 (-3) 

1 . 32C-6) 

4 

6 

58.7 

1 .97177 (-1 ) 

6.46C-5) 

3  =  4,  M1  = 

M2  =  5*  M3 

=  =  15,  Tr  =  446.5 

sec . 

4 

7 

71  .5 

4.67654(-4) 

7  - 1 4 (-8 ) 

1 

6 

60.5 

9 .82626 (-3) 

1 . 25 (-6) 

0 

6 

60.3 

1 . 98789 (-1  ) 

6 .46 (-5 ) 

3  = 

■  6,  M1  -  M2 

=  m3  =  m4  = 

M5  “  M6  "  5>  TR  "  305 

.9  sec. 

[8 

6 

48.6 

4 .59736 (-4) 

7.76(-8) 

•  3 

6 

43.6 

9 . 64348 ( —  3 ) 

1 . 58 ( - 6 ) 

3 

6 

45.9 

1 .93787 (-1 ) 

6 .73( -5) 
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